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Chapter 1 


INTRODUCTION 

1.1 General Background 

The concept of an Earth and Ocean Model, in itself, is 
nothing new. Every development program involving missiles, 
satellites, geodesy, etc., must take into account one or more 
geophysical influences, each of which must be mathematically 
modeled. The degrees to which these influences are considered 
and related within the model are variable but important 
factors. 

Most models ancillary to development programs have been 
restricted in scope and generally confined to one-time appli- 
cation. Traditionally, this is a result of limitations in 
accuracy and scope of measuring systems and consequent inabil- 
ities to isolate interrelated geophysical phenomena. It has 
long been known in theory that separated points on the Earth's 
surface were in relative oscillatory motion, but these motions 
were in the decimeter range and, from the viewpoint of measur- 
ing potential, could not be discriminated from the "solid" 
Earth. 

The advent of artificial satellites and of sophistication 
in electromagnetic measuring devices and transportable time 
standards afforded the potential and imminent realization of 
intercontinental geodetic measurements to precisions in the 
centimeter range. This was formally recognized and documented 
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at a seminar conducted in August 1969. The report (Kaula 
1970) ) defined the existing and anticipated State-of-the-Art, 
and was a basis for an Earth and Ocean Physics Application 
Program (EOPAP) . 

The EOPAP as described in the contract has been devel- 
oped to study and utilize the discipline of Earth dynamics 
for the benefit of science and mankind. Earth dynamics 
include solid-Earth and ocean dynamics, which are concerned 
with physical motions, and distortions of the Earth and phys- 
ical state of the ocean. Solid-Earth dynamics include forces 
responsible for earthquakes, tidal waves, volcanic eruptions, 
mineral differentiation and mountain building. Ocean dynamics 
embraces ocean circulation and physical state of the ocean 
surface. A thorough understanding of Earth dynamics is funda- 
mental to intelligent management of the Earth and its impact 
on our scientific knowledge and application. Progress in 
solving the problems of environmental management and allevia- 
tion of natural disasters will be quite limited until a better 
understanding has been achieved of physical mechanisms that 
respond to these forces. Earth and ocean dynamics are enor- 
mously complex, involving important interactions between 
components in different parts of the globe, and in the atmo- 
sphere, oceans and solid Earth. This precludes any probabil- 
ity of finding solutions in terms of simple theoretical des- 
criptions. Determination of the probability of catastrophic 
natural events will almost certainly have to be based on 
strongly empirical numerical computer models that require, 
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as operational inputs, very large numbers of current synoptic 
data. This computer model was intended to identify the need 
and mission requirements for flight systems to interface with 
the user's requirements; to interrelate various geophysical 
and oceanographical phenomena, and to systematically catalog 
and store data for ready recall and application. The tasks 
of developing and applying a satellite system for a World- 
wide network of data points will be impossible without such 
a method of correlating and analyzing the collected informa- 
tion. 

1.2 Objectives and Scope 

The major objective of this study was to identify and 
formulate an Earth and Ocean Model (EOM) for use in the EOPAP 
and to develop a system that would supply information to the 
user for applying physical data gathered on a World-wide 
basis. It was intended that this would lead to the develop- 
ment of theoretical and empirical formulae which could be used 
to determine physical forces and influences of geophysical and 
oceanographical phenomena. 

The EOM is effectively a modular structured system of 
computer programs utilizing Earth and ocean dynamical data 
keyed to finitely defined parameters. The model is an assem- 
blage of mathematical algorithms with an inherent capability 
of maturation with progressive improvements in observational 
data frequencies, accuracies and scopes. The programs provide 
facilities for accepting a broad spectrum of data, not only 
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for evaluation of equation coefficients but for incidental 
storage and retrieval for later program extensions. To 
insure a simple yet practical method for updating geophys- 
ical constants and parameters, the values which were selected 
are listed in Appendix A and, for the computations, are iso- 
lated in a block data subroutine. This provides a ready means 
for positive, one-point entry for updating program parameters. 
The main program interfaces data input with appropriate module 
subprograms and controls interrelated routines. 

It was originally intended that the EOM would include 
the following geophysical modules: 

a. Plate tectonics; 

b. Earth rotation and polar motions; 

c. Ocean currents; 

d. Ocean tides; 

e. Earth tides; 

f. Magnetic fields; 

g. Gravity; 

h. Atmospheric motion; 

i. Sea state; however, during the investigation 
it was established that ocean currents and sea states were not 
amenable to practical modeling - the former, since there 
•appeared to be considerable disagreement in the literature 
(and between oceanographers) of critical transport phenomena, 
and the latter, because of a real-time dependency on short- 
term variations in the atmospheric mechanism. Required 
computing facilities would be far beyond the scope of the EOM 
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program. ^ Atmospheric motion was restricted to its effect on 

Earth rotation and polar motion.^ 

A geodetic system where relative coordinates of widely 

separated points on the non-rigid, deformable Earth's surface 

are instantaneously determined to centimeter accuracies 

requires not only a recoverable and stable reference system 

but also a means for identifying the short- as well as long- 

period variations from the equilibrium Earth figure. The 

requirement for the former was the principal concern at an 

International colloquium on reference coordinate systems. ^ 

The following keynote was presented by Lundquist (1974) : 

"The current need for more precisely defined reference 

coordinate systems arises for geodynamics because the Earth 

can certainly not be treated as a rigid body when measurement 

uncertainties reach the few centimeter scale or its angular 

equivalent. At least two coordinate systems seem to be 

required. The first is a system defined in space relative to 

appropriate astronomical objects. This system should approx- 

^ The matter of ocean currents and sea states was discussed 
with the Director, Space Sciences Lab., MSEC, and the COR. 

It was decided to delete the ocean current and sea state 
modules. 

^ Atmospheric motion is determined by processing continuously 
available meteorological data. Global limitations (political 
and physical) preclude proper distribution and frequency of 
observation stations. This radically limits the effective- 
ness of atmospheric motion determinations for other than 
localized treatment, and are consequently unsuitable for the 
EOM application; however, the synoptic effects of atmospheric 
motion are reflected in Earth rotation variation and polar 
motion and can be evaluated from statistical analysis of these 
effects. (Approved at meeting described in above.) 

International Colloquium on Reference Coordinate Systems 
for Earth Dynamics, Torun, Poland, August 1974. 
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imate an inertial reference frame, or be accurately related 
to such a reference, because only such a coordinate system is 
suitable for ultimately expressing the dynamical equations of 
motion for the Earth. The second required coordinate system 
must be associated with the non-rigid Earth in some well- 
defined way such that the rotational motions of the whole 
Earth are meaningfully represented by the transformation pa- 
rameters relating the Earth system to the space-inertial 
system. The Earth system should be defined so that the dynam- 
ical equations for relative motions of the various internal 
mechanical components of the Earth and accurate measurements 
of these motions are conveniently expressed in this system. " 

It was accordingly decided to include a precession 
module as an extension of the polar motion module and to 
treat the Earth rotation in a separate module. The modules 
which were actually modeled and their program names are as 
follows: 


a. 

Plate 

tectonics (EOMPTA) 

b. 

Earth 

rotation (EOMERA) 

c. 

Polar 

wobble (EOMPWA) 

d. 

Precession (EOMPMA) 

e. 

Ocean 

tides (EOMOTA) 

f . 

Earth 

tides (EOMETA) 

g- 

Magnetic fields fEOMMFA) 

h. 

Gravity (EOMGMA) 
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The rather large excursions and somewhat periodic secular 
motion of the Chandler polar path was investigated in some 
detail. It was demonstrated that both could be attributed to 
internal Earth's mass redistributions marked by large earth- 
quakes. 

1.3 Outline of Report 

Geophysical considerations relating to development of 
the EOM modules are described in Chapter 2; assumptions, 
simplifications, and acceptance of specific theories are 
discussed in sufficient detail to define and validate adopted 
approaches and associated provisions. Chapter 3 contains 
mathematical and physical details involved with assembling 
the individual EOM modules, as well as limitations and poten- 
tial for further development. Chapter 4 presents UNIVAC 1108 
system FORTRAN V programming details for the EOM modules to 
include a comprehensive indexing of computer materials. 
Applications of the EOM and factors bearing on future exten- 
sion are discussed in Chapter 5, including the investigation 
of the power source for the polar wobble mechanism. Chapter 
6 provides a summary of the EOM program and appropriate recom- 
mendations and conclusions. 



Chapter 2 


GEOPHYSICAL CONSIDERATIONS AND CONCEPTS 
2 . 1 General 

The Earth's interior has been and still is a largely 
undetermined geophysical variable; however, seismic tech- 
niques (Takeuchi, et al. (1967)) have provided bases for 
assumptive but fairly conclusive determinations of the inter- 
nal structure. Supporting the atmosphere and hydrosphere is 
a relatively thin crust enclosing, first, a mantle, and then 
a core. Crust and mantle are solid while the core is taken 
to be liquid (outer core) with the deeper part (inner core) 
generally considered to be solid. The mantle-crust boundary 
is marked by a significant change in constituent material 
called the Mohorovicic discontinuity after its original dis- 
coverer. The crust thickness varies over the Earth from about 
5 kilometers under oceans to about 70 kilometers beneath 
mountains. The average continental crust thickness is about 
35 kilometers. There are several published models for the 
internal properties of the Earth, but Earth Model B (Bullen 
(1963)), shown in Table 2.1, is representative. The outer 
core is assumed to act as a self-exciting dynamo. The dynamo 
theory, first suggested by J. Larmor in 1919, was worked out 
in detail by W.M. Elsasser and E.C. Bullard (Rikitake (1966)) 
following World War II. The principal portion of the geomag- 
netic field is due to this influence of the core and the 
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Table 2 . 1 

Properties of Earth 

Model B 


Depth 

P 

P 

k 

y 

g 

33 

3.32 

0.009 

1.16 

0.63 

985 

80 

3.36 

0.025 

1.22 

0.66 

986 

80 

3.87 

0.025 

(1.40) 

(0.76) 

986 

200 

3.94 

0.071 

(1.58) 

(0.83) 

985 

400 

4.06 

0.150 

(1.92) 

(0.99) 

983 

600 

4.18 

0.231 

(2.61) 

(1.34) 

980 

1000 

4.41 

0.400 

3.37 

1.78 

976 

1400 

4.63 

0.58 

3.96 

2.03 

976 

1800 

4.84 

0.76 

4.58 

2.25 

982 

2200 

5.03 

0.96 

5.26 

2.48 

997 

2600 

5.22 

1.16 

5.89 

2.71 

1010 

2700 

5.27 

1.22 

6.13 

2.81 

1042 

2898 

5.57 

1.33 

6.40 

2.97 

1069 

2898 

9.74 

1.33 

6.4 

0.0 

1069 

3500 

10.60 

1.95 

8.2 

0.0 

937 

4000 

11.16 

2.42 

10.1 

0.0 

815 

4500 

11.63 

2.85 

11.6 

0.0 

647 

4982 

12.00 

3.22 

12.1 

0.0 

607 

5121 

15.4 

3.33 

13.6 

(3.2) 

573 

6371 

18.1 

3.95 

16.4 

(5.0) 

0 

Depth in 

kilometers; 

p is the 

density in 

g/cm^ ; p is 

the 


pressure xlO’^ dynes/cm^ ; k is the incompressibility and y the 
rigidity, respectively, in the same units as p; g is the 
acceleration of gravity in cm/sec^. Discontinuities indicated 
at 80, 2898, and ca. 5050 kilometers represent crust-mantle- 
outer core-inner core. 
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remainder, minute and highly localized, due to remanent mag- 
netizations in the crust. There is strong evidence that the 
core rotates at a different rate than the mantle-crust system. 
This is manifested in a westward drift of the geomagnetic 
field and in variations in Earth rotation rate which may be 
due to momentum conservation for varying core-mantle relative 
angular velocities. 

There is, of course, much more known about the Earth's 
surface. A continuing mystery, which seems to have been 
fairly resolved only in the last decade, concerned paleoecolo- 
gists' efforts to account for inconsistencies in distribution 
of land masses and of ancient flora and fauna. The first 
scientifically reasonable explanation, in terms of continental 
drift, was given in 1915 by A. Wegener (Wegener (1924)).** His 
theory was not immediately accepted mainly for lack of a sci- 
entifically plausible explanation for the driving forces 
required to move continental masses. Continental drift theo- 
rists were eventually opposed by polar wandering theorists,^ 
with a third group accepting a combination of the two (Hapgood 
(1970)). There is now a consensus for continental drift in its 
scientifically structured form known as plate tectonics. 

Plate tectonics are concerned with the tectonic activities 
of an Earth's surface divided into a mosaic of rigid, shifting 

English translation of the original "Die Enstehung der' 
Kontinente and Ozeane" , Brunswich, Germany (1915). 

® Polar wandering assumes the outer shell to be intermit- 
tently shifting about the poles moving continents collec- 
tively. Continental drift treats continental motions as 
individual . 
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plates. There are 14 major plates whose boundaries are asso- 
ciated with a seismically active variety of characteristic 
features such as rift valleys, oceanic ridges, mountain belts, 
volcanic chains and deep oceanic trenches. Oceanic ridges 
mark regions where adjacent plates are being pulled apart and 
new surface material is formed from cooling magma. Oceanic 
trenches are formed by one of the two plates plunging steeply 
into the Earth below the other. The total effects of these 
actions are to maintain a constant Earth's surface area and 
to displace surface features in varying degrees with respect 
to an absolute geodetic datum. 

It has been estimated from seismic data that plates are 
from 70 to 150 kilometers thick and, accordingly, not related 
to the mantle-crust differentiation marked as the Mohorovicic 
discontinuity. Best evidence indicates that they form a rigid 
layer overlying a weaker, hotter layer which becomes increas- 
ingly viscous with depth. The supporting surface is called 
the asthenosphere . Plate boundaries move at relative veloci- 
ties of from 1 to 10 centimeters per year. Since this motion 
is on the surface of a sphere the relative motions of any two 
plates can be described by a rotation around a single axis 
passing through the sphere's center. Given the absolute 
motion of any one plate with respect to a. set of fixed points 
within the mantle, the absolute motion of all the plates can 
be determined from knowledge of the individual relative plate 
motions. Appropriate processes for this are described by 
Dewey (1972) and Solomon and Sleep (1974) . 
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The sea-level surface of the Earth very nearly conforms 
to an ellipsoid of revolution with equal equatorial axes. If 
the axis of rotation is along a principal axis, the angular 
momentum vector is coincident. It was suggested by Euler 
(1765) that if a rigid figure did not rotate about a principal 
axis, a free nutation would result as a function of the prin- 
cipal moments of inertia. For Earth values this would have a 
period of about ten months. It was reported by Chandler 
(1891) , after an extensive study of latitude variation period- 
icity, that a ten-month period did not exist; instead, there 
was an annual term and a second term of about fourteen months. 
It was shown by Newcomb (1892) that this unexpected increase 
from ten to fourteen months was due to mobility of the oceans 
and an elastic yielding of the Earth. This non-annual nutata- 
tion or polar wobble® is now known as the Chandler wobble. 

The angular momentum vector, in the absence of external 
forces on the Earth would retain a fixed orientation in iner- 
tial space. It lies within the plane formed by the axes of 
rotation and figure, very nearly coinciding with the former. 
The rotating Earth is subject to gravitational forces of Sun, 
Moon and planets, and the subsequent gyroscopic effect results 
in precessional motions’ of the angular momentum vector. The 

® Polar wobble is considered to be the trace, of the rota- 
tion axis with respect to some axis on the Earth's surface 
to which latitude stations are referenced. Polar wobble 
is thus manifest in terms of latitude changes. 

Precession is a measure of the motion of the angular momen- 
tum vector through the celestial sphere. The axis remains 
inclined at approximately the same angle to the ecliptic 
plane and is manifest in declination changes for the stars. 
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process is complicated by non-coincidence of axes of rotation 
and figure, deformation of a yielding Earth by gravitational 
forces, and seismically induced mass redistributions in the 
heterogeneous structure of crust and mantle. 

The Earth's rotation rate is influenced by external 
torques, changes in inertial moments, wind stress and core- 
mantle coupling. External torques include effects of Earth 
and ocean long-period tides, gravitational torques on the 
equatorial bulge, and the solar wind. Changes in inertial 
moments include seismically induced internal mass redistribu- 
tions, surface level changes brought about by polar ice and 
sea level fluctuations, and atmospheric mass transport. Var- 
iations in the rotation rate are monitored by astronomic time 
observations and, more precisely, by utilizing such techniques 
as the Very Long Base Inter ferometeric (VLBI) procedure for 
observing on extragalactic radio stars (Gold (1967)). 

The entire Earth's surface is in continuous, relative, 
periodic and systematic motion, and there are no direct means 
for selecting and maintaining an absolute geodetic reference 
system which can be readily and positively recovered. Since 
the reference system must be available at the surface, the 
best solution is to find an approach which minimizes undetect- 
able displacements. An ideal system would have the origin at 
the Earth's center of mass and three or more axes aligned to- 
ward selected extragallactic radio sources. The origin is 
readily recoverable through appropriate satellite observa- 
tions. The same would be true of the inertial axes using 
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VLBI except for the lack of a sufficiently accurate means for 
relating them to the stellar system. The present choice is 
to set the origin at the Earth's center of mass, the z axis 
toward a point known as the Conventional International Origin 
(CIO)®, the X axis orthogonal to the z and pointing toward a 
prime meridian defined by the Bureau Internationale de I'Heure 
(BIH) ® in Paris. The y axis completes the orthogonal set, the 
positive direction pointing towa:rd 90° E. longitude. It 
should be noted that geophysists utilize this right-handed 
system, while astronomers choose the positive y direction 
toward 90° W. longitude. 

Rochester (1973) has tabulated the spectrum of changes 
in Earth's rotation and associated geophysical mechanisms. 
These are presented in Tables 2.2 and 2.3, respectively. 
Elements in columns A, B and C of Table 2.2 are discussed in 
detail in Sections 2.10, 2.14 and 2.13, respectively. 

The Earth's elastic properties can be effectively sum- 
marized in a set of dimensionless parameters known as Love 
numbers. Love (1909) introduced the numbers h and k to char- 
acterize the various aspects of Earth tides. His theory was 
keyed to the concept that since the tidal potential can be 
adequately represented by a second-degree spherical harmonic 
function, all deformations due to this potential may be eval- 

® The CIO is defined as the average position of the true 
celestial pole from 1900 to 1905 as determined from adopted 
latitudes of the International Latitude Service (ILS) lati- 
tude stations. 

The BIH determines changes in adopted observatory longi- 
tudes necessary to maintain an absolute prime meridian that 
is not necessarily the same as the Greenwich meridian. 
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Table 2.3 Mechanisms with Effects now Distinguishable 

on the Earth's Rotation 


Mechanism 

Effect* 

Sun 

Gravitational torque 
Solar wind torque 

A, B7, Cl, C3c 
C2c (?) 

Moon 

Gravitational torque 

A, B7, Cl, C3d 

Mantle 

Elasticity 
Earthquakes 
Solid friction 
Viscosity 

Bl, B3-4, Cl-2a, C3c-d 
Bl, B3 
B3 (?) , Cl 
C2a 

Liquid core 

Inertial coupling 
Topographic coupling 
Electromagnetic coupling 

A2-3, B2, B6 
C2b-c (?) 

A4 (?) , B3, C2 

Solid inner core 

Inertial coupling 

B2(?) 

Oceans 

Loading and inertia 
Friction 

Bl, B3, B5, C2a 
B3 (?) , Cl 

Groundwater , 

Loading and inertia 

B4 

Atmosphere 

Loading and inertia 
Wind stress 
Atmospheric tide 

B4 

C2c, C3a-c 
Cl 


* 

Letters and numbers refer to Table 2.2 
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uated by multiplying the harmonic function by the appropriate 
Love number; e.g., the Earth surface is raised by ^-h, where U 
is the surface potential and g is the acceleration of gravity. 
The consequent gravitational potential is given by U*k. The 
horizontal analogue of the Love numbers, was suggested by 
Shida and Matsuyama (1912) . It is; generally treated as an 
additional Love number. Application of the Love numbers are 
discussed in detail in Section 2.8. 

Although the point has been made that the Earth is not a 
rigid body and must be treated apart from rigid body theory, 
the departures from rigidity are sufficiently small to be 
considered as perturbations to the rigid body solutions. 

Since these solutions are well known, a combination of pertur- 
bation technique and Love number application provides a com- 
pletely adequate treatment. An effective approach was pre- 
sented by Munk and MacDonald (1960) and is adopted in this 
study. 

The theoretical influences of lunar and solar potential 
on an elastic Earth are conveniently expressed in terms of 
spherical harmonics. Laplace (1796-1825) separated the poten- 
tial into three groupings; he called these the tides of the 
first, second and third type which are uniquely representable 
as zonal, tesseral and sectorial harmonics, respectively (See 
Section 2.15) . (These correspond to long-period, diurnal and 
semi-diurnal periods, respectively.) The total tide repre- 
sents the sum contributions of all its components at a speci- 
fied time and place. Component arguments can be expressed in 
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terms of six independent astronomical variables. There are 
31 principal components (Melchior (1966)) with relative am- 
plitudes ranging from 910 to 3. The practical majority of 
the tide can be determined by using three long-period, three 
diurnal and three semi-diurnal components. 

A reasonable solution for the Earth tide can be obtained 
using second-degree spherical harmonics if Love numbers are 
utilized to account for elasticity, and some lag factor is 
introduced to account for delay in Earth-potential response 
time. 

The situation is much more involved for ocean tides 
which include the open ocean (pelagic) tides as well as the 
more familiar ones in the vicinity of land masses. The indi- 
cated approach for the Earth tide is essentially a static 
one which can not be applied to the ocean tide. The principal 
problems are Earth's rotation (Coriolis effect) and the limit- 
ing effect of ocean depth on tidal wave speed. Ignoring the 
Coriolis forces for the moment, a direct response of the ocean 
to lunar or solar gravitational attractions requires that the 
bulge travels at the same speed as the attracting body. If it 
does so, the result is known as the direct tide; however, the 
tidal wave speed is dependent upon water depth and the bulge 
can fall behind, leading to a 90 degree phase shift and an 
inverted tide. Since the apparent speed of the attracting 
body is a function of latitude this would, by itself, lead to 
radical discontinuities in ocean levels and mass transport. 
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This, of course, does not occur, but the matter of dynami- 
cally evaluating the pelagic tide is a problem v/ith no pos- 
sible direct solution. Assumptions must be made of average 
water depths, optimized boundary configurations, transport 
processes, etcetera. Boundary conditions are the observed 
tides at or near land masses. In practice, solutions are 
made for separate components, each of which is manifest as 
an amphidromic system (See Section 2.16). 

2.2 Reference Frames and Transformations 
2.2.1 Geographic Frame; 

The first considered reference frame is a set of axes 
rotating with the Earth. There are several possibilities, 
but the most feasible and readily implemented system is a set 
of orthogonal axes defined by the coordinates of a network of 
observatories. In this sense the axes are fixed to the solid 
Earth and rotate with it. The existence of a deforming Earth 
and its effect on a rigid set of axes imposes a requirement 
that the observatory coordinate definition be modified to in- 
clude provisions for variations due to oscillatory and secular 
motions. 

The z axis, which very nearly represents the rotation 
axis, is defined as the Conventional International Origin 
(CIO) . The CIO was adopted in 1967 by the International 
Astronomical Union (lAU) and the International Union of Geod- 
esy and Geophysics (lUGG) . It set the pole to the average 
position defined by the five stations of the International 
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Latitude Service (ILS) for the period 1900-1905. The Bureau 
Internationale de I'Heure (BIH) converted their system to the 
CIO by adjusting the BIH polar paths to those of the ILS for 
the years 1964-1966. The x axis is orthogonal to the z axis 
and in the plane of a zero meridian determined by longitude 
(time) observations of the BIH. The y axis is orthogonal to 
the xz plane and considered here as positive in the direction 
of 90° E. longitude. 

The BIH utilizes data from its latitude and longitude 
stations and from the Dahlgren Polar Monitoring Service 
(DPMS) to provide continuous observatory updates for long- 
period and secular Earth motions. Provisions for higher fre- 
quency motions are made in the data reduction processes of the 
individual observatories. 

The coordinate origin is the Earth's center of mass. 

This is readily recoverable through satellite observations. 

2.2.2 Inertial Frame: 

The inertial reference system is made up of three or- 
thogonal axes oriented with respect to the ecliptic and ver- 
nal equinox of some selected epoch. Here the epoch of choice 
is 1950.0 or the Julian date, 2433282.5. The inertial refer- 
ence system is related to the Earth-fixed rotating system 
through the three Euler angles. The relationship of the two 
frames is shown in Figure 2.1. The positive directions of the 
Euler angles 0, (j) and are indicated in the figure. The dots 
over the Euler elements indicate the time derivatives. The 




/ Vernal E< 
(1950 


Equator of CI< 


BIH definition of 
Greenwich meridia: 


sumption that the CIO axis remains 
actice, the CIO axis is replaced b 
and the CIO equator is defined by 
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inertial axes are X,Y,Z; rotating axes are x,y,z. 

The kinematical relationships of the Euler angle deriv- 
atives are developed by adding details of Figure 2.2 to those 
of Figure 2.1. The axes x,y,z of Figure 2.2 repeat those of 
Figure 2.1. The additional axes are the angular momentum axis 
OH and the instantaneous rotation axis OP. (Note that H lies 
on the great circle arc zP.) It can be seen from the two 
figures that 


0 ) = 0 ) sin Y cos T , 

1 

0 ) = 0 ) sin Y sin T , (1) 

2 

(jj = 0 ) cos Y / 

3 


where is the polar rotation vector, the components 1,2,3 
representing x,y,z, respectively; o) is the rotation vector 
magnitude and y and F are angles indicated in the figure. 

The rotations o)^ which determine the motion of the CIO axis 
with respect to the instantaneous rotation axis represent 
the same motion as the Euler angle derivatives. These rela- 
tions are found by resolving each of the vectors into 

its components along the CIO axis. The result is 


U) 

1 

= 

^-sinSsincf) 

-COS(j) 

0 

w 

2 

= 

-sin6cos(f) 

sine}) 

0 

U) 


COS0 

0 

1 

3 







The inverse relationship is 




= 

”-csc0sinc() 

-CSC0COS(}) 

O' 


(U 

1 

* 

e 

= 

-cos<p 

sin(j) 

0 


0) 

2 


= 

cot0sin(() 

COt0COS(}) 

1 


to 

3 


Note that since 0 remains on the order of 23°. 5 no mathemati- 
cal difficulties are presented by the first of Eq. (3). 

The transformations between the system and the 
system obey the following relations: 


X 


i 



a 





(4) 

(5) 


where 
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a 


ij 


cos())cosii;-sin(j)sirnjjcos0 cos(J)sirnjj+sin(j)COSi|;cos6 
-sin(t)cosi|j-cos(|)sin4)Cos6 -sin(j)sini{j+cos(t>cosijjcos6 
-sini|isin0 cosi|jsin0 


-sin(J)sin0 

-coscj)sin0 

COS0 


( 6 ) 


2.2.3 Geodetic and Associated Frames: 

A geocentric system of coordinates can be expressed by 
the equations. 


X 


(v+H) cos4) cosX 

y 

= 

(v+H) cos(() sinX 

z 


_[v(l-e^)+Hj sin(j) 


where: 4) is the geodetic latitude; \ is the longitude (positive 
in the eastern hemisphere) ; H is the height of the point above 
the ellipsoid, and 


V = 


e‘ = 


(l-e^sin^cf)) ^ 
a^-b^ 


.2 


( 8 ) 

(9) 



(10) 


with a and b the semi-major and semi-minor ellipsoid axes, 
respectively. The inverse of (7) is 


X = tan-^ (^) 

X 


( 11 ) 


and the iterative solution of 
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<|) = sin 


- 1 


lv(l-e^)+H 


( 12 ) 


2 

where 0 = sin" (— ) is introduced into v for the initial 
determination. The 4) in v must be successively updated until 
the increment in the determined <}) becomes negligible. When 
H=0 the exact solution is 


<J) = tan 


- 1 


z sinX 


y(l-e^) 


= tan 


- 1 


z cosX 
x(l-e") 


(13) 


The geocentric latitude 4/ is given by 


4< = tan“^ ^(1-e^) tancl)j , 


(14) 


and for the radius vector p. 


cos(j) ^ h 
‘cos'l'cos ( 4>-fp) 


= a( = a(l-f sin^cj)) 


(15) 


where the flattening f is 


f = 


a-b 


(16) 


for points on the ellipsoid surface. For points at elevations H, 

Pg = Pg + H . (17) 

Given geocentric rather than geodetic coordinates, the geo- 
centric coordinates are 


X ' 


p„ cosip cosX 
n 

Y 

~ 

p„ cosip sinX 

n 

z 


pjj sini|; 
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The accuracy of (18) is improved by substituting ip„ for 

rl 

where ij^„ is given by 

n 


'H 


= tan 


- 1 


v (l-e^)+H 

v+H 


tanc}) 


(19) 


The reduced latitude 3 for a point on the ellipsoid is given by 


3 


tan" 


b . . 

— tand) 
a 



( 20 ) 


Ellipsoidal coordinates u,0,X are related to rectangular co- 
ordinates by (where 0 = 90° - 3# the reduced latitude) 


X 


(u^+E^ ) ^cos3cosX 


a cos3cosX 

Y 

= 

(u^+E ^) ^cosBsinA 

= 

a cos3sinX 

z. 


u sin3 


b sin3 


(21) 


where E is the linear eccentricity, 

E = (a^-b^)^ . 


( 22 ) 


Figure 2.3 Latitude Definitions 
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Altitude corrections of latitudes are given by (with H in 
meters) , 

A(j) = 1”.71 H sin2(j) x 10““ ^ AS - (23) 

The latitudes are defined in Figure 2.3. In Figure 2.3A, 
the radius vector p is the length OP. The line nP is a 
normal to the ellipsoid surface. OE is the equatorial line 
intersected by the plane ONP. In Figure 2.3B, the line Ob 
equals (u^+E^)^, the line Oa equals u, and these two radii 
represent the semi-major and -minor axes of the ellipsoid 
surface. P, in both figures, is the observation point. 

2.2.4 Astronomic Frames : 

The astronomic coordinate systems of interest here are 
the ecliptic and right ascension systems for some specified 
epoch or instant. The ecliptic system locates the astronomic 
body in respect to the ecliptic and vernal equinox; the right 
ascension system, in respect to the Earth’s equator and vernal 
equinox. Ecliptic coordinates are latitude, 6, and longitude, 
A. Right ascension coordinates are declination, 6, and right 
ascension, a. The two systems are related by the equations 


cos6cosa = cos3cosA , (24) 
sin6 = cos3sinAsine + sin3cose , (25) 
cos6sina = cos3sinAcose - sin3sine , (26) 
cos3cosA = cosScosa , (27) 
sin3 = -cos(Ssinasine + sinScosE , (28) 
cos3sinA = cos6sinacose + sin6sine , (29) 
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where e is the obliquity of the ecliptic. Given the ecliptic 
parameters the value of 6 is obtained from (25); a is 

obtained after dividing (26) by (24). Given the right ascen- 
sion parameters 6,a,e, the value of 3 is obtained from (28); 
the value of X is obtained after dividing (29) by (27) . 


2.3 Vector Relationships Applicable to Plate Motions 

Rotational components of an ellipsoid normal with coor- 
dinates (j),X, are given by 


(i) 


COS(j) cosX 

X 


0) 

= 0) 

cos(j) sinX 

y 




sin(f) 


(30) 


u) = is the rotation magnitude and o)^ the rotation vector. 

The displacement of a point moving on a plate rotating on an 
axis with coordinates (p,X, at a rate can be established 
from the diagrams in Figure 2.4. The relationships associated 
with the diagrams are given with the figure. From these. 


i. e. , 


’Ax’ 


3 

3 

1 

0 

1 


X 



z y 



Ay 

= t 

X 

3 

1 

0 

N 

3 


y 

Az 


0 

3 

3 


z 



_ y X J 






X t 


f 


(31) 


(31' ) 


where is the completely antisymmetric permutation tensor 

13 K 

(refer to text on tensor calculus; e. g. , Hotine (1969)). 
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From Eq. (7) with H=0 and setting (l-e^)=l (noting that 
(l-eM=0.9933) , 


- 

X 


cos4) cosX 

y 

= V 

cos(J) sinX 

z 


sin4) 


(32) 


dX = 


d(j) = 


X = tan * (^) 

X 


4) = sin ^ (^) 


3X , ^ 9X , 

55 57 


dz 

vcoscj) ' 


} 


(33) 


x^+y^ 


(xdy - ydx) , (34) 

(35) 


where v is considered constant. Introducing (31) into 

(34) and (35) (noting that lim = lim = — ) , 

z-^0 v(l-e^) 

zt 


dX = 0) t - 

z 


(XO) + yw ) , 

2 • 2 X y 

X +y ^ 


d(j) = — (yo) - XO) ) t 

^ z X y 


(36) 


(37) 


In Eqs. (36) and (37) , t represents the interval of time in 
years from the epoch (time at which positions were assumed to 
be fixed) to the time of observation. The quantities are 
obtained from adopted values of the absolute plate motions 
which are prescribed here in Appendix A. 

Note that — = tan<j)secX, — = tan4)cscX, cosX = r- , 

^ ^ (x^+y^)^ 

^ — T- . ■ From Eq. (36) and these, 

(x"+y=^)^ 


sinX 
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AA = 


0) 


W (1) 

- 2 cosX sinA (— + — ^) 


CO - tand) (co cos A + co sinA) 
z X y . 


( 38 ) 


In a like manner Eq. (37) can be written as 


Ad) = (co sinA - co cosA) t 

^ X y 


(39) 


2.4 Complex Parameters for Excitation Functions 

The various geophysical effects upon the Earth's motion 
will be treated in terms of excitation functions. They are most 
simply treated in complex form. All forms of 'i' in the fol- 
lowing are complex: 


'V = 

cosat + 

sinat = 


(40) 

4 / = 

4 ,+j ^i(at+A‘^) 

+ If”! e’ 

(at-A ) 

(41) 


+ 




T = 

(arg H* ) 

= 

(arg ) 

(42) 

4 /^ = 

4 /+ + 4 ,- 

4-^ = 

i(4,f _ 4 ,-) 

(43) 

II 

+ 

jg(4,c _ j^4,Sj 

4»“ = 

+ i4'®) 

(44) 

f = 

■ ■ ^ — — (ellipticity 

4-'^ + 4 )" 

or flattening) 

(45) 


The major axis lies in east longitude Js(A^+A ) and has a mag- 
nitude I +j’l' I . It is occupied at such a time that 
at = % ( A -A^) . 
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2 . 5 Interpolation Formula 

Interpolation of non-linear tabular data (e.g., ephem- 
eris data) is by Gauss's forward interpolation formula 
(refer to Hildebrand (1956) for the general form) where the 
desired interpolated value is given by 


f ^ +^{A +^^{A + 

o+nw o 1 2 2 3 3 


+il^{A }}}}} 


where: A = f^-fo 

A — 1 

2 

A = f 2~3f i+3f Q-f-i 

3 

A = f 2~4f 1 +6f p-4f_ j+f_2 

4 

A = f 3-5f 2+lOf j-lOf o+5f_i-f_2 

5 


(46) 


and the subscripts are obtained from the Julian dates. 


J.D. 


Eph. 


J.D. 


Obs. 


Subscript . 


2.6 Spherical Harmonics 

The Legendre functions for the spherical harmonics used 
in the EOM may be any one of four types: the regular form 

denoted by P or (See Heiskanen and Moritz (1967)); an 
alternative form denoted by p^ or p^ (See Jeffreys and Jef- 
freys (1962)); a partially normalized form denoted by P^ or 
^ (See Schmidt (1964)), and a fully normalized form denoted 
by P' or P'’'” (See Heiskanen and Moritz 1967)). These dif- 
ferent forms are related as follows: 
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Pn = 

P 

n 

(47) 

m 

Pn = 

(n-m) ! m 
n 1 *^n 

(48) 

^n = 

P 

n 

(49) 

^ = 
n 

,2 (n-m) I , *5 
' (n+m) 1 ^ n 

(50) 

P' 

n 

= (2n+l)^ P^ 

(51) 

p.n' 

n 

= 

2(2n+l^ ^ P"' 

(52) 


The Legendre expression can be calculated using the equa- 
tion (modified from Munk and MacDonald (I960)), 


^m 

n 


= 2“%! (n+m) ! (1-t^)"'^^ 5 


n-m 


■ (t-l)^'^~^(t+l)^ 

(m+r) ! (n-r) I (n-m-r) Ir! 


where t=sincl). 


2 . 7 The Normal Potential ^ ® 

The normal potential can be expressed using only zonal 
harmonics because of the rotational symmetry. The potential 
series is of the form (where p=sin(})) , 


V = “ + A 
r 


P (y) 


+ A 


2 j, 3 


4 ^5 


A2n = (-1) 


n kME 


2n 


2n+l 


(1 - 


2n me ' 
2n+3 3q 


(54) 

(55) 


E = (a^-b^)* = ea = e'b , 


(56) 


q^ = 2 (^ e' 3 e'5 + e*’ 


■3-5 


5*7 


7*9 


) . ( 57 ) 


Source is Heiskanen and Moritz (1967). 


(53) 
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Eq. (55) can also be written as 

A = (-1)" 3GME^" 

^2n ^ ^ (2n+l) (2n+3) 


(1 - n + 5n — ) 
ME^ 


(58) 


If the potential V is written in the form. 


V 


[l - f: ' 

n=l 


r 


(59) 


n. 


then the J (noting = -Qla J^) are given by 


/ i\^tl *5 2n 

'^2n " 


(2n+l) (2n+3) ^ ^ ^2 ' 


From the general form of Eq. (54) , 

00 n 


V = i: ^ 

n=0 m^ 


^ R^(0,X) „ S^(0,X) 

An -~n+T ®n " Iri+I ' 


dM 






(60) 


(61) 


> (62) 


The solid harmonics and are simply homogeneous 

polynomials in x, y and z. For example. 


r^S^ = 3r^sin^0sin2X = 6r^ sin^OsinXcosX = 6(r sinOcosX) (r sinOsinX) , 
2 

= 6xy . 
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The moments of inertia with respect to the x, y, z axes are 


A 

B 

C 




)(y'"+z'M 

dM 

|(z ' ^+x' ^ ) 

dM 

J(x'"+y'2) 

dM 

[ X ' y ' dM 



I (63) 


From Eqs. (62) and (63) , setting A=B, 


A = GM 
0 


A = A^ = 
1 1 


B^ = A^ = 
1 2 


B^ = A^ = 0, 
2 2 


A^ = G(A - C) , bI= 


(64) 


(Note that GM for the Earth is here taken as 3.98603 xlO^** MKS, 
The gravity formula in closed form is 


Y = 


GM 

a(a^sin^6 + b^cos^3)^ 


e'q' e 'q' 

,.,m Ov . 2 n , /I It' 

(1+2 — - — )sin^B+(l-m-g — — 


) cos^ 3 


At the equator (3=0) , 


e'q' 

Y = •SMd-m-E ; 

'a ab 6 q 


( 66 ) 


At the poles (3= i90°) , 


GM in ® ’ Q 

Y, = —(1+^ 


(67) 


Note that q^ is given by Eq. (57). For 


' 6(^ 


— © ' ^ ^ — e ' ** + — — — e ' ^ — 

5 5*7 7-9 ® 


) , ( 68 ) 


m = 


a)^a^b 

GM 


) 


. (65) 
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With f=— , f*= - — - 
a Y, 


IS 


, the rigorous form of Clairaut's formula 


f + f* = 


(jj^b 


(1 + 


e ’q’ 


2 q, 


(69) 


The gravity formula (65) can be written in series form as 


Y = Y. 


1+ (-f+|m+jf ^-^fm+^^m^ ) sin^(})+ (-^f ^+^fm) sin**(^j , ( 70 ) 


or 


Y = Y, (1 + f ,sin^4)+f . sin*‘(|)) , 


(70 ' ) 


with 

fj = -f+^+jf ^-^fm+^pm^ , 

f = -if^+^fm 

4 2 2 

Eq. (70') can be written as 

Y = Yg (1+f 4sin^2(J)) , 


^(71) 


(70") 


where f* is known as the gravity flattening. 

The normal gravity at elevation H above the ellipsoid is 

Y = Y+— • • • 


= Y 


l--(l+f+m-2fsin24»)H+-HM . 

Si o 


(72) 
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2.8 Love Numbers and Applications 

Love niimbers may, in theory, be taken to any degree n. 

The measure of application is the nature of the disturbing 
potential. Since the Earth is effectively subject only to a 
second-degree potential. Love numbers of consideration are h, 
k and H . The matter has been studied by a number of research- 
ers,^^ but the most detailed investigation appears to be that 
of Takeuchi (1950) . Love number determinations are made from 
observations of seismic phenomena, and tidal induced tilts, 
variations of gravity and position. The k coefficient can 
also be determined from observations of polar wobble period. 
Estimates are subject to contributions of Earth's core, mantle 
and yielding oceans. Since l(ove numbers are indicators of 
elastic characteristics, values should vary randomly for dif- 
ferent locations; in addition, a latitude dependency has been 
suggested by Kaula (1969) based on analysis of published data. 
There is, however, no strong reason to systematically charac- 
terize variabilities in Love numbers, and, for the purposes of 
the EOM, values of h=0.59, k=0.29 and £.=0.07 have been adopted. 
The Chandler period is given by Kaula (1968) as 

T = Atka.^aV(3G) ^ 

C-A-kw^a V(3G) 

'' These include Jeffreys (1970), Kaula (1968), Kozai (1965), 
Lamb (1917) , Love (1911) , Longman (1966) , Melchior (1966) , 

Munk and MacDonald (1960) , Nishimuri (1950) , Shida and Mat- 
suyama (1912) , Takeuchi (1966) , Vicente (1961) . 
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The products of inertia due to rotational deformation (Munk 
and MacDonald (I960)) are 


^i 


^ = I6 . ^ - ^^6.^) 

1 3G 1 3 1 


(74) 


where I=4(C^+C^+C^ ) is the Earth's inertia, a is the mean 
radius and G is the Earth's mass. Note that the effect of 
elasticity in (73) is to lengthen the period and in (74) to 
induce products of inertia. 

The luni-solar potential is effectively given (See 
Section 2.9) by 

W = - ^ —(1 - 3 cosH) . (75) 


For a non-deformable Earth, the horizontal components of the 
luni-solar force are 


f = — 

^ Q 1 -» 


, aw 

1 2 


R 30 


9W, 


f. = 


R sin0 ax 


^(76) 


where 0 is the colatitude and X the east longitude. The cor- 
responding vertical deflections are 



^(77) 


where is the meridian and n ' the prime vertical components. 
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respectively. The spatial displacements for an elastic Earth 
are 


V 


w 




g 36 


gsinS 3X 


^ (78) 


in the meridian and prime vertical directions, respectively. 
The total disturbing potential is (l+k)W 2 and a static oceanic 
tide should have an effective height (l+k)W 2 /g; a tide gage, 
however, indicates variations between ocean tide and crustal 
motion, the latter being given by ^ 2 * ^ observed tidal am- 

plitude will, therefore, be the difference between these two, 
or 


^2 h ^2 

(l+h)-f - = (i+k-h)— = Y— 

g g g g 


(79) 


Given as the initial Earth's potential, gravitational 
acceleration is 


g 


o 


GM 


8V 

o 

3R 


But the Earth's potential, deformed by the luni-solar action, 
is 


V 


W + W ' + 

2 2 


3V 

o 

“5r" 


+ V 

o 




and the variation of g will be the derivative with respect to R 
of the additional potential. 
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- g 


= 


’ + 


3V 

o 

“5F 


(80) 


Here ' represents the additional potential due to the 
Earth's deformation (in the form S^/R^) while 



Performing the indicated differentiations in (80) , 

W , - 3W 8W2 

Ag = 2^(l+h-|k) = (l+h-|k)^ = . (82) 

The angle between the normal and the deformed surface 
is of .interest. The equation for the latter is 

p = R + h — . (83) 

The flattening of this surface can be calculated by relating 
its equation with that of an ellipsoid of revolution with mean 
radius a and ellipticity e; i.e., 

p' = R[l+^(|-cos2e)] , (84) 


where 9 is the colatitude. From (84) and the value of g for 
a homogeneous sphere of radius R and mean density d, or 

4 

g=jirGdR, the deformed eccentricity is given by 


2 _ 


.GhRM 


gr 


(85) 


The angle a of the normal to this ellipsoid with the unperturbed 
ellipsoidal normal is 



a 


sin20 


p M 2 


3 GMhR 

^ gr^ 


sin29 


1 h 
g R 


3W2 

TF • 


This quantity added to the pendulum deviation, 
gives the observed effect as 


Rg 3 9 


f 


1+k-h 

Rg 


3W^ 


3W 

= ^ ^ 
Rg 3 6 


( 86 ) 


The crust deformation is in the same direction as the force 
and thus is partly compensating. The factor y is the ratio 
of the observed deflection amplitude to the calculated theo- 
retical static amplitude. 

Deflections cause a variance in astronomic coordinates 
The deflection is not related to the crust but to the direc- 
tion of the Earth's rotation axis in inertial space. Eqs. 
(77) and (78) express the respective components of vertical 
deflection and surface deformation. The deformation will al 
cause a supplementary disturbing potential kW^ which must be 
added to in the expressions. Finally, the deflection of 
the vertical with respect to the Earth's axis will be 


? = (l+k-£) 


1 9^2 

Rg 39 


Rg 30 


3W, 


n = il+k-l) 


Rgsin0 3 A 


A 

RgsinO 3X 


^( 87 ) 



2.9 Tidal Force 


The tidal force acting on a unit mass at point P due to 
the attraction of an extraterrestrial mass M is 


F. = GM 


x.-X. X. 

11 1 




where M is the mass of the attracting body. Expanding the 
denominator of the first term and assuming jx^| = R << |x^| = r, 

F. = -—(X. — -^i) = -—(X. - — - OS^x . ) . (89) 

1 J.3 1 ^2 J.3 1 r 1 


Eq. (89) can be resolved into its horizontal and vertical 
components which are 


F(h) = |gm— sin2c , 


F(v) = GM— (l-3cos^<;) 

^ 3 


( 91 ) 
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It can be seen that Eqs. (90) and (91) can be derived from 
the potential. 


. 7 GM 2 r 

(3 cos^C - 1) , 

^ r 3 


by taking the derivatives, 

1 

R 3c 


and 


F(v) = - 


3W 

2 

Tr“ 


(92) 


\ 


2.10 Precession and Nutation 
2.10.1 Louiville Equation: 

The Eulerian equations of motion in a coordinate system 
rotating with angular velocity relative to coordinates 
fixed in inertial space are 


N. = H. + 
1 1 


e . . , 
xjk 


(93) 


and are torque and angular momentum vectors, respectively, 

and e. is the completely unsymmetric permutation tensor. 

1]K 

Following Munk and MacDonald (1960) , we separate the angular 
momentum into two parts: 


H. 

1 


. (t)w^ + h^(t) 


(94) 


where 


^i 




p(x,x^5? - x.x^)dV 

K X X 


(95) 
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is the inertial tensor, p is the density and 6^ is the Kron- 
ecker delta. The second right-hand term of ( 94 ) denotes a 
relative angular momentum. 


'i = ■3'' ' 


( 96 ) 


due to the velocity u^ relative to the system. Introducing 
( 95 ) and ( 96 ) into ( 93 ) we obtain 


N. = ^ (C. .0)^ + h. ) + e. .,0)^ (C^^Wo + h^) 
1 dt ij 1 ijk £ 


( 97 ) 


the Louiville equation. 

2 . 10.2 Dynamic Equations for the Precession: 

In terms of rectangular coordinates x^' of M' (the 
external attracting mass) in the CIO coordinate system (with 
BIH prime meridian) , the couple, N^, exerted on the Earth is 


N. = e..,x^'F^ = - e..,x^'V^U 
1 13k 13k 


( 98 ) 


Introducing these into ( 97 ) , 


C. .(1)^ + C. + h. + e. .. 0)^ (C^^a)„+h^) = e. ..x^'V^U 
13 13 1 13k £ 13k 


Assume, for the moment, that the relative angular momenta, h^ 
are included in the angular momenta, and these are constant 
Then ( 99 ) becomes 


C. .0)^ + e. = 

13 13k £ 


- e . x^ ' V^U 
13k 


( 100 ) 


( 99 ) 


We next make the temporary assumption that there are no products 
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of inertia. Then setting j=A, C^^=B, C^^=C, ^^^=0, i^j 
and expanding (100) , 


• C-B 

+ -^2^3 




• 3U V 


to, - to 

2 B 13 


1, , au 

B ^^1 Sx^ 


3^^ 


( 101 ) 


• B-A 
to. + to to 

'"3 ^ ^^12 


1/ , 9U 


, 3U V 


The right-hand terms of (101) can be written in terms of 
Euler angles and we obtain 



, C-B 1 

A 2 3 A 

sin (|) , a 

3U 

3(J) 

3U. 

- COS(j) 

3U 

36 

^2 

C-A 1 

- — 7— <0 0 ) = — 

B 1 3 B 

sin9 

3U 

3(J) 

3U) 

+ sin4) 

3U 

36 


(102) 


. ^ B-A 1 3U 

“3 + -c^i “2 = c 3* 


We can set A=B, whence 9U/3(j)=0 because of the symmetry, and 


tOj + 


“2 - 


C-A 

_ sincj) 

^ _ 

COS(p 

3U 

A 2 3 

A sin9 


A 

36 

C-A 

_ COS({) 

^ + 

sine}) 

3U 

_ UJ UJ 

A 1 3 

A sin6 

dip 

A 

36 


(103) 


(Oj = constant. 

2.10.3 Precession and Nutation for a Rigid Earth: 

The complexity of the problem indicates that the best 
solution would be to treat the precession using a rigid Earth 
and making such corrections as are necessary for the elastic 
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yielding in a separate process. This provides a highly accurate 
representation and will afford firm bases for correcting for 
departures from rigidity. We follow to some extent the work 
of Woolard (1953) . 

The first assumption is that there are no products of 
inertia; secondly, that the moments about the two minor axes of 
figure are equal. Eq. (103) is thus applicable. Differenti- 
ating the first two of Eq. (3) , we eventually obtain 


0 = - 


Coo jSin6 


dtp 


Co), dt 


(4»sin9) + 


Coo 


-l(l0COS0 


(104) 


and 


ip = 


CoOjSinO 


3U 

30 


Coo 3 sine 


0 + 


^ ijj^cose 


Cu) 


(105) 


Eqs. (104) and (105) can be utilized in an approximate 
form by considering only their first right-hand terms. This 
gives a close approximation of the axis of figure. Small ad- 
ditions to account for the effect of the unused terms are 
appreciable but can be easily obtained by an iterative solution, 
the convergence being quite rapid. 

The force function U for the gravitational action exerted 
on the Earth by a external rigid body with mass M' and principal 
moments of inertia A',B',C, is 


U = 



^ 2r^ 2r^ 


where and are moments of inertia of Earth and M' , respec- 
tively, about the line r joining the C.M.s of the two bodies. 
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In the case of luni-solar forces, effects of extended terms 

and of the A',B',C inequalities are so small we can take 

S'+B'+C'=3I and 
r 

U = GM' (- + A+B±^^) . (106) 

"" 2r^ 

Setting the Earth's coordinate axes along the principal iner- 
tial axes and denoting the coordinates of the C.M. of mass of 
M' , in this system, by we have with x^x^=r^, 

I = — (Ax' 2 +Bx' ^+C x' . (107) 

r j,2 1 2 3 

Since r does not depend upon the orientation of the Earth's 
principal axes, it is independent of the Euler angles with 
respect to which U is differentiated in the equations of mo- 
tion; therefore, all terms in U which depend only upon r and 
constants can be disregarded. Setting A=B in (106) , we 
obtain 

U = - 3GM'— xl^ . (108) 

2r^ 

In terms of the system, the coordinates of the disturbing 
body are 

X = r cosB cosA 

1 o o 

X = r cosB sinA 

2 o o 

X = r sing , 

3 *^o ' 


W(109) 
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where and are the latitude and longitude, respectively, 
of M' referred to the fixed mean equinox and ecliptic of 
specified epoch. From Eqs. (4) and (109) , 

x' = r(sin9cosB sin(X -iji) + cosOsinB ) • (110) 

3 o o o 

Introducing (110) into (108), differentiating, and introducing 
into (104) and (105), as appropriate, (utilizing only the 
first right-hand terms of these latter two) we obtain 

9 = -3GM* {sin9cosB sin(A -ij)) +cos0sinB }cosB cos(A -i|;) (111) 

r^Co ,3 ° ° ° ° o ^ ^ ^ 

and 

>p = -3GM* sinOcosB sin(A -ij^) +cos9sinB }{cot9cos6 sin(A -i|j) 

r^Coo 3 o o o o o 

-sinB } . (112) 

o 

The coordinates B »A , of the Sun and Moon, referred to 

o o 

a fixed mean equinox and ecliptic of some epoch are obtained 
from solar and lunar theories; but these theories give the 
coordinate expressions referred to a moving mean equinox and 
ecliptic of date, and it is therefore necessary to transform 
Eqs. (Ill) and (112) to these coordinates. Through an appro- 
priate process (refer to Woolard (1953)), the coordinate re- 
lationships between B^»A^ and B/A (the latter being the coor- 
dinates of date) are 
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cos3qCos ( = cosB{cos ( X-A+ n j -ij;) +2sin ( A-A) sin (n ^ -ijj) sin^hiT j 

+sin3siniTjSin(nj-ij;) , 

cos^^sin ( = cos6{sin(X-A+n^-<ji) -2sin (X-A) cos sin^SsTTj } 

-sinBsinTT jcos (IT , 

sin3Q = cos3sin ( X-A) sim; j+sin3cosn j , 

where the elements are defined in Figure 2.6. 

An alternate and much more simple procedure may be 
used in developing the quantities 3 q/X^ from 3(t), X(t) for 
Eqs. (Ill) and (112). We have, for an epoch of 1S50.0, 

T = (J-D. - 2433282. 5)/36525 

e = 23°. 445787 , 

o 

e = 23°. 445787 - 0°.0130125T - 0°. 00000164T^ 

+ 0°.00000050T^ + Ae , 

where Ae is the nutation in the obliquity; then, 

6 = sin” Mcos3sinXsine + sin3cose) , 
a = cos” * (cos3cosX/cos6) , 

= z -0".791T^ - 0”.0013T^ 

0 = - 2004". 255T + 0".426T^ + 0".0416T^ 

q = sin6{tan6 + cos (a+C^) tanJsG } , 

Aa - y = tan”^{qsin (a+C^) / l-qcos(a+CQ) } r 
= a + Aa , 

6^ = 6 + 2tan”^{tanis6seci5 (Aa-y) cos [(a+?^) +*5 (Aa-y) ] } , 

See Am. Ephem. and Naut. Aim., NHO, 1976, pp. 536 and 552. 
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z = - 2304". 948T - 0".302T^ - 0".0179T^ 

U = z + Co » 

3 = sin~*(-cos6 sina sine + sin6 cose ) , (114) 

o o o o o o 

= cos”^ (cos6^cosa^sec3o) - (115) 

A less rigorous transformation from 3 (t) , X (t) to 3o» 
X^ is obtained from the expressions, 

3q == 3 + b sin(X+c) , (116) 

X^ = X + a - bcos ( X+c) tan3 , (117) 

where 


a = -1°.396319T - 0°.000308T^ , 

b = -0°.013076T + 0°.000010T^ 
c = 5°. 59258 - 1°.15473T - 0°.00015 t2 

Setting 

K = - 3GM' ■ 

r'Co)3 

P = sin6cos3 sin(X -tlj) + cos0sin3 
o o o 

Q = cos3 cos(X -\1^) , 

o o 

R = cos3 sin(X -il;)cot0 - sin3 , 

•^o ' o ^ o 


>( 118 ) 
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Eqs. (Ill) and (112) can be written as 

0 = KPQ , dll') 

= KPR . (112 ' ) 


Noting that 


ijj = K(PR + PR) , 

P = KP{Qcos0sin (X^-i{;) 
- Qsin0sing^ , 

Q = KPRcosB sin(X -\p) 
o o 

R = - KP{Qsin(X -ip) + 
o 


(119) 

- Rsin0cos(X -tJj) }cos3 1 

o o 


/ 


W(120) 


*5Rsin26cos ( X^-iJ;) }cos3^csc^ 9 . 


Introducing Eqs. (Ill'), (112'), (119) and (120) into Eqs. (104) 

and (105) , we obtain the bases for iterative solutions of (104) 
and (105) , 

0 = KPQ + (PR + PR)sin0 + 2KQRP^COS0} , (121) 

CO)3 


= KPR - 


AK 


Coo 3 sin 


(PQ + PQ - JsKP‘'R‘'sin20) , (122) 


where the underbars are introduced to discriminate from the 
provisional values 9, ip, used in forming the right-hand terms, 
which will be successively updated in the iteration process. 

Final values of 0 (t) and i|^(t) are obtained by integrat- 

• • 

ing updated ^(t) and ij^(t) quantities. The remaining Euler 
angle, 4> r is obtained from the expression. 
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= Tp + T , (123) 

where <=y»^ is the ephemeris value for the Julian date of the 
midnight preceeding the first observation date, and t is 
given by 

T = 360°. 98565 (J.D. , . - J.D.^) 

(obs) o 

2.10.4 Correction for a Deformable Earth: 

The effects of the luni-solar force on a deformable 
Earth are considered in Section 2.15. Values which substan- 
tially compensate for the rigid Earth assumption made in this 
section are derived. 

2.11 Plate Tectonics 

Plate tectonics theory has been extensively developed 
over a relatively short period, with the investigative pro- 
cess being ultimately directed toward discovering a scientif- 
ically valid mechanism for the plate motions. Relative plate 
motions, manifested as continental drifts, are fully accepted 
phenomena and have been (for at least the ten largest plates) 
fairly well evaluated. The matter of absolute motions; 
i.e., motion with respect to the fixed portion of the mantle, 
say, the asthenosphere , presents a more 'illusory problem, 
since there are no proven means for referencing to a fixed 

Refer to Chase (1972), LePichon (1968), Minster et al. 
(1974) , Morgan (1968) . 




triad in the inner mantle. A number of authors^** have pre- 
sented theories for absolutely relating one or more plate 
motions to fixed interior references. The consensus links 
absolute motion determination with the plate driving mech- 
anism, and most writers have tried to bring them together 
as a total system. 

Plate boundaries are fairly well delineated by system- 
atically tracing seismic activities. Maps, such as those 
described by Barazangi and Dorman (1969) provide indications 
of plate demarcations and directions of motion. Kaula (1975)^^ 
outlines the boundaries of thirteen plates in a system of 
great circle arcs. These plates and their boundaries were 
incorporated in the EOM. A fourteenth, the Caribbean, was 
added, with boundary based on Malfait and Dinkleman (1972) . 

A number of absolute plate motion systems were evalu- 
ated for application to the EOM. The system of choice was 
one of the models proposed by Solomon and Sleep (1974). The 
selected Model A fixed the absolute motion of one plate based 
on a general concept that there was a uniform drag beneath 
all plates. This simple model agreed remarkably well with 
more sophisticated ones (Morgan (1973), Minster et al . (1974)), 

whose velocities were calculated from the hypothesis of fixed 

^ These include Burke and Wilson (1972), Clague and Jarrard 
(1973), Duncan et al. (1972), Kaula (1975), Mckenzie et al . 
(1974), Minster et al. (1974), Sleep (1975), Solomon and Sleep 
(1974), Solomon et al. (1975), Turcotte (1975), Wilson (1965). 

Referred to by the author in a computing program. Program 
details and boundaries are on computing tape. (Refer to 
Section 4.10.) 
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hot spots, and there is no basis for a more involved system, 
with the present state of knowledge.^® The absolute motions 
of the remaining plates were calculated from the relative 
motion data given by Solomon and Sleep (1974) . 

Plate parameters are the coordinates of the rotation 
pole, and the rotation rate or, more simply, only the compon- 
ents of the rotation vector. Plate names, identifying numbers 
and rotation parameiters are given in Table 2.4. An example 


Table 2.4 Tectonic Plate Parameters 


0 )^ (10“’deg/yr) 


Number 

Plate Name 

Abbreviation 

0) 

X 

0) 

Y 

w 

z 

1 

African 

AFRC 

1.3483 

-2.4585 

2.8415 

2 

Antarctic 

ANTA 

-0.2660 

-0.4475 

2.6859 

3 

Arabian 

ARAB 

4.8787 

-1.0630 

4.5355 

4 

Cocos 

COCO 

-7.3499 

-15.2714 

7.5398 

5 

Eurasian 

EURA 

-1.2936 

-1.1125 

-1.4502 

6 

Indian 

INDI 

5.4483 

2.0823 

4.4787 

7 

Nazca 

NAZC 

-2.2031 

-4.6138 

6.0178 

8 

North American 

NOAM 

0.1281 

-1.7780 

-0.3030 

9 

Pacific 

PCFC 

-1.2867 

2.8497 

-6.4968 

10 

South American 

SOAM 

0.1281 

-1.7780 

-0.3030 

11 

Bering 

BERI 

0.1281 

-1.7780 

-0.3030 

12 

Philippine 

PHLP 

3.2305 

1.0976 

-6.0047 

13 

Somali 

SMLI 

1.3483 

-2.4585 

2.8415 

14 

Caribbean 

CARI 

-0.4474 

-0.2825 

0.5532 


* ® Any model can be readily introduced into the EOM by 
making the desired changes in the rotation parameters. 
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for updating plate parameters follows: 

Given the absolute rotation pole parameters for the 
Pacific plate (PCFC) as (|)=-64°.3, A=114°.3, w=7 . 21x10“ Meg/hr ; 
given the pole parameters of the African plate (AFRC) relative 
to PCFC (assumed fixed in motion) as (j>=57*^.6, X=-63*^.6, o)=11.06 
xlO” ’deg/hr; find the absolute rotational components of AFRC: 
From Eg. (30) , 


Absolute 

PCFC: 

0) = 
X 

-1.2867, 

03 = 2.8497, 

y 

03 = 

Z 

-6.4968. 

Relative 

AFRC: 

0) = 

X 

2.6350, 

03 = —5 . 3082 , 

y 

03 = 

Z 

9.3383. 

*^i (PCFC) 

■^‘^i (AFRC) ■ 

03 = 

X 

1.3483 

03 = -2.4585 

y 

03 = 

Z 

2.8415. 


The inverse procedure (rotation components to pole pa- 
rameters) is given by^® 


0) 


1 h 

0) . (i) 

1 


(positive root, only) , 


0) 

<t> = sin“ ^ (-^) , 

(0 

(0 

X = tan“ ^ (— ^) 

CO 


^(124). 


The general boundaries of the plates are shown in Fig- 
ure 2.7. A more exact plot is given on the map attachment to 
the computing materials (See Section 4.2.3). in most cases, 
identification of the appropriate plate for a station of given 
coordinates will be quite simple; however, when a station lies 


Note that, for all poles 
gitudes are negative. When 
1X1000 . 


, latitudes are positive and lon- 
(0 is negative, |X|> 90 O; otherwise 


/ 
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near a plate boundary, the boundary description may be too 
general, and a definite check must be made to identify the 
pi^te . For this reason, and because of the imprac— 
ticability of mathematically delineating the highly irregu~ 
13-1^ ^nd often changing boundaries, it was decided to reguire 
a hand— ident i f ication " process for all points; thus avoiding 
problems for critically located sites. 


Figure 2.7 Tectonic Plate Boundaries 


(Adapted from Solomon and Sleep (1974)) 
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2.12 Atmospherics Applied to Polar Wobble 

The polar wobble exhibits an annual component due to 
atmospheric variations. Specific influences can be calculated 
from appropriate pressure field measurements, if these are made 
at sites uniformally distributed over the Earth; unfortunately, 
this appears to be physically and politically impossible, and 
there are extensive regions from which needed data can not be 
obtained i The impact of inadequate coverage was noted by 
Wilson (1975) : he compared atmospheric excitation functions 
published by a Soviet investigator, with functions derived by 
others and concluded that the rather significant differences 
were due to restricted availability of Siberian meteorology. 

Atmospherics, like any phenomena conforming to a spheri- 
cal surface, can be mathematically treated, most conveniently, 
using spherical harmonics. This has been done by a number of 
authors. The process is complicated, however, by the need to 
maintain orthogonality of the functions. Kaula (1967) de- 
scribes statistical analyses of data distributed over a sphere, 
and rigorous requirements to insure validity. Data gaps cause 
non-orthogonality and require analytical measures involving 
covariances; mere substitution of null values, followed by 
least-squares determinations or other devices, may result in 
inordinately specious harmonic coefficients. Although harmonic 

Includes Awade et al. (1975), Eliasen (1958), Eliasen and 
Machenhauer (1965) , Ellsaesser (1966) , Graham (1955) , Hassan 
(1960) , Haurwitz (1940) , Moses (1974) , Neatan (1946) , Platz- 
man (1960) . 
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methods can be advantageously employed for limited area cov- 
erage, their applications to world-wide atmospherics do not 
appear to be feasible due to the limited scope of our present 
meteorological observing system. 

Straightforward methods, for analyzing pressure field 
motions, probably provide a more satisfactory means for deter- 
mining atmospheric influences. Seasonal variations in air 
mass distribution, which contribute to the annual component 
of polar wobble, can be estimated from the ground level pres- 
sure field. Seasonal wind variations do contribute to changes 
in the earth rotation rate, but, as shown by Munk and MacDonald 
(1960) , these appear to have little or no effect on polar wob- 
ble. The data gap problem still exists, and there appears to 
be little likelihood of approaching a 50 percent representation. 
Since the actual effects can be determined rather positively 
from spectral analysis of long-period latitude observations (See 
Section 2.14), atmospheric analysis appears to offer little more 
than academic interest in the polar wobble application. 

2.13 Earth Rotation 

Earth rotation rate, referred to the inertial frame, 
affords a general measure of sidereal time. The fundamental 
unit of sidereal time is the mean sidereal day, the interval 
between two successive transits of the mean vernal equinox over 
a given meridian, corrected for polar motion and certain peri- 
odic irregularities in earth rotation rate. The more readily 
utilized solar time is based on the local hour angle of a 
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fictitious sun. Solar time, when referenced to a defined 
prime meridian, is termed universal time (UT) . There are 
four UT systems in use: UTO is observed UT; UTl is UTO cor- 

rected for the BIH definition of the prime meridian (corrected 
for polar wobble) ; UT2 is UTl corrected for periodic variations 
in earth rotation; UTC (coordinated universal time) is the 
system used in radio time signal transmission, and is an 
approximation to UTl, obtained from international atomic time^® 
(TAI) by the relationship, UTC = TAI - 15®. 

Rotation rate is continuously monitored by observatories 
of the BIH, which issues a monthly report of five-day raw and 
smoothed time and polar wobble data; these include values for 
UT2-UTC and UTl-UTC. More accurate determinations will become 
available with implementation of new techniques, most impor- 
tantly, that of VLBI. 

Principal variations in earth rotation rate are caused 
by tides, polar wobble and zonal wind circulation. Lesser, 
periodic influences, are, due to atmospheric pressure and hydro- 
logical variables (snow, ground water, vegetation, etc.). 

Secular changes in rotation rate are due to core-mantle coupling 
and effect of lunar torque on the tidal bulge of the Earth. 

Tidal effects include fortnightly, monthly, semi-annual 
and annual terms; polar wobble influences the rate over a four- 
teen-month period; zonal winds add semi-annual, annual and, 
intermittently, near-biennial contributions. Atmosphere pres- 

TAI is directly based on the data of atomic clocks. 




sure and hydrological variables occur at semi-annual and annual 
frequencies. The core-mantle coupling effect is manifest in a 
secular variation of the geomagnetic field (See Sction 1.17); 
it is a consequence of angular momentum conservation during 
variations in rotation rates between core and mantle. The ef- 
fect of a lag in developement of the tidal bulge (caused by 
imperfect elasticity and, thus, energy dissipation) is a torque 
by the Moon on the Earth, slowing the earth rotation. 

Tidal effects can be theoretically calculated, but the 
meteorological measuring situation noted in Section 2.12 applie 
and estimates of zonal wind influence are subject to consider- 
able uncertainty; evaluations of atmospheric pressure and hydro 
logical effects are highly conjectural; secular contributions 
are fairly observable. 

Influences on earth rotation can be characterized by 
dimensionless excitation functions (See Sections 2.4 and 2.14). 
Lambeck and Cazenave (1973) give a number of these functions 
which are listed in Table 2.5. 

Polar wobble does not have a direct effect on earth 
rotation. Astronomic coordinates of the time observatories 
must be corrected for the effect of polar wobble to insure 
proper orientation with the CIO and BIH prime meridian. The 
fortnightly and monthly tidal terms do not lend themselves 
well to the format used for the excitation functions of the 
other terms; the lunar motion introduces a complication. 
Excitation functions for these (Munk and MacDonald (I960)) are 

See, also, lijima and Okazaki (1972) . 
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Table 2.5 Excitation Functions 


Period 
T (Months) 

Coefficients 


A 

B 

Geophysical Factor 

6 

-0.051 

0.260 

Zonal wind (altitude 0-30km) 

6 

0.050 

0.010 

Zonal wind (altitude 30-60km) 

6 

-0.0006 

0.0000 

Atmospheric pressure 

6 

0.003 

-0.005 

Ground water 

6 

0.161 

-0.060 

Earth tides 

12 

-0.344 

-0.131 

Zonal wind (altitude 0-30km) 

12 

-0.024 

-0.010 

Zonal wind (altitude 30-60km) 

12 

-0.0032 

0.0027 

Atmospheric pressure 

12 

0.017 

0.022 

Ground water 

12 

-0.028 

-0.001 

Earth tides 

24 

-0.077 

0.029 

Zonal wind (altitude 0-30km) 

24 

-0.001 

-0.003 

Zonal wind (altitude 30-60km) 

6 

0.162 

0.205 

Total for 6 month period 

12 

-0.382 

-0.117 

Total for 12 month period 

24 

-0.078 

0.026 

Total for 24 month period 

Excitation function. 

'I's = A 

cos — ^ + B sin where t is the 

number of 

months from Jan 0 , 

1958, and coefficients A and B are 

given in 

units of 10 

- 8 



- = 0 . 16xl0“®cos (2s-N) 

+ 0.38xl0“®cos2s (125) 

and 

~'^Mm ~ 0 . 20xl0“®cos (s-p) , (126) 

where Eq. (125) is the excitation function for the fortnightly 
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term and Eq. (126) , for the monthly term. The arguments are 

s = 270°. 434358 + 13°. 1763965268d - 0°.001133T ^ 

o o 

+ 0°.0000019T^^ , (127) 


N = 259°. 183275 - 0° . 0529539222d + 0°.002078T ^ 

o o 

+ 0°.000002T ^ 
o 


(128) 


p = 334°. 329653 + 0° . 1114040803d - 0°.010325T ‘ 

o o 

- 0°.000012T ^ 
o 


(129) 


where s, N and p are the mean longitude of the moon, longitude 
of the mean ascending node of the lunar orbit, and mean longi- 
tude of the lunar perigee, respectively. and d^ are 


defined as 


T = d 736525 
o o' 


(130) 


and 


d = J.D. , . . - 2415020.0 

o (obs) 


(131) 


The total excitation function (Table 2.5 data and Eqs. 
(125) and (126)) is (in units of 10“®) 

27Tt 2TTt 27Tt 

= 0.162 cos^ +0.205 sin^ - 0.382 cos^ 

3 D b X ^ 

2tt 1" 27r1" 2Trt 

- 0.117 sin^ - 0.078 cos^ + 0.026 sin^ 

- 0.16 cos(2s-N) - 0.38 cos2s - 0.20 cos(s-p) . (132) 


The equation of motion is (See Eq. (140) 

t 


m = -p; (-o)c - h + N dt) = 'I' 

3 0)C 3 3 3 3 3 


IS 


( 133 ) 
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where the right-hand side is represented by the arbitrary 
evaluation in Eq. (132) . (Elements of (133) are defined in 
Section 2.14.) From Eq. (140), 

, (134) 

“3 

where, from Eq. (137) , m, = — - 1, is the "direction cosine - 1" 
relating the rotation axis to the axis of figure. (See Fig- 
ure 2.2 and footnote to Figure 2.1.) 

The effect of core-mantle coupling is more difficult 
to characterize. Munk and MacDonald (1960) discuss the 
marked westward drift of the geomagnetic field and suggest 
the possibility of comparing observed westward drift with 
mantle rotation. A number of additional authors^® have in- 
vestigated the apparent correlation between geomagnetic and 
astronomical time observations. The consensus is that the 
correlation is positive and valid, it is difficult to numeri- 
cally evaluate, and the matter bears further investigation. 
This item is discussed in more detail in Section 2.17.^^ 


These include Aoki (1969), Bondi and Littleton (1948), 
(1953), Bullard et al. (1950), Elsasser (1950), lijima and 
Okazaki (1972), Kahle et al. (1967), Malin and Saunders 
(1973) , Munk and Revelle (1952) , Nagata (1965) , Rochester 
(1960) , (1968) , (1973) , Runcorn (1954), Siran (1969), Steen- 
beck and Helmis (1975) , Triet ^974) , Vestine (1952) , (1953) , 
Vestine and Kahle (1968) , Yukatake (1962) , (1972) . 

The effect of the core-mantle coupling on polar wobble, 
as well, has been reported by several investigators; this 
will be discussed here. 
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Lambeck and Cazenave (1974) discuss uncertainties in 
several of their suggested excitation function values (shown 
in Table 2.5.) These estimated errors, coupled with greater 
ones for other functions, indicate that an earth rotation 
model could be subject to considerable error; however, the 
literature^ ^ presents evidence that uncertainties are sig- 
nificantly smaller than magnitudes. In addition, the EOM 
design affords a ready input facility as more valid data be- 
comes available. The matter of correcting astronomical longi 
tude observations for polar wobble effects (discussed above) 
is covered in Section 2.14, and provisions are made for such 
corrections in the computing modules. 

A final note concerns the possible effect of solar 
winds on the earth rotation. Coleman (1974) suggested that 
the solar wind could measurably affect the rotation rate; 
Hirshberg (1972) and Hines (1974) point out that the effects 
would not be significant; Gribbin and Plagemann (1973) claim 
a variation of 10 msec (jump in AT-UT2) between August 7 and 
8, 1972 during a great solar storm, confirming their predic- 
tions of such results; O'Hora and Penny (1973) saw no such 
result and, since a storm of the observed magnitude provided 
such minimal effects, feel that there are good grounds for 
believing that changes in the length of the day are induced 

^^Refer to Brouwer (1952), Challinor (1971), Fliegel and 
Hawkins (1967), Frostman et al. (1967), Guinot (1970), 
lijima et al. (1964), lijima and Okazaki (1972), Lambeck 
(1975) , Markowitz (1970) , Mintz and Munk (1951) , (1954) , 

Munk and MacDonald (1960) , Munk and Revelle (1952) , 

Rochester (1973) , Sidorenkov (1968) , (1973) , Woolard (1959) . 
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by some other mechanism. Evaluation of these reports indicates 
little likelihood for variations in earth rotation due to so- 
lar wind or flares. 


2.14 Polar Wobble 

For convenience, the mathematical development for the 
polar wobble and Earth's rotation are discussed in this section. 
We follow Munk and MacDonald (1960) in the general approach. 

The basic equations for the polar wobble and rotation are the 
Liouville equations given in (97) , which, in component form, are 


♦‘i* 31c 21c 

+ h^ + a)j^(o)2C “‘^3^ ) + ^*^2^3 " ^3^2 ' 

«2 = =2j“^ * =2j“^ * i'2 ^ + “ 3^1 - “l"2 , 

•T •*!• 21c lie 

+ C^jW-* + h^ + o)j^(Wj^C ) + w^h2 - ^^2^1 


H135) 


Consider the perturbation scheme where^ in the moment of inertia 
tensor. 



A+Cii 

*^12 

•=13 ] 

c. . = 

13 

^21 

^■^^22 

=23 


^31 

^32 

C+C 33 


(136) 


A, B and C are the moments of inertia referred to the principal 
axes and c^j are small corrections and products of inertia for 
i=j and respectively. Let 
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= wm^, = wnij,, (»3 = oid+nij). 


(137) 


1 i I ^ 

where w=|a)j^a) j . Introducing (136) and (137) xnto (135) and 
setting B=A, 


m, = 


1 2 


= 


(D^ (C-A) 
1 


2 .2 


0)^ (C-A) 
1 


(“"Ci3 + C0C23 

+ wh^ 

+ hj 

+ A(joiTi2 

- Nj) 

((d^C22 “ 

+ (Dh2 


- Aionij^ 

+ S 3 ) 

33 ^3 " ^3^ 






X(138) 


c . . h. 

It can be seen that m. and are small, dimensionless 

C 1 (OC 

quantities whose products and squares can be neglected. We set 


^ - C~A 

r A 


(139) 


Introducing (139) into (138) , noting the vanishing products, 

^2 1 - ^2 ^2 

- — + m, = ir::— (wc,, + c^^ + h, — )= - 4), , 

1 AO^ 13 23 1 (Jj U) r 

^1 3 ' 

— + m„ = — ((DC-- - c,_ + h- - — + — )= (})- y.(140) 

2 13 2 0) (D 2 , 


m. 


= - 1(C + 1 f 

C'^33 (D (dJ 


N^dt) = 


The right hand terms represent the dimensionless excitation 
functions, <j)^, and are determined by geophysical observations; 

/ 

the left-hand terms are determined by astronomical observa- 
tions . 

The first two of (140) represent the polar wobble equa- 
tions of motion. They are more conveniently handled in complex 
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form, with 


m = m^^ + im 2 <t) = 4)^^ + i<t )2 t 


and written as the single equation. 


• O' . X 

1 — + m = d) 

a 

r 


(141) 


The third of (140) is the equation of motion for the earth 
rotation. Quantities N^, h^ and in (140) are qualified 

as follows: 

a. The initial time t marks the instant the rotating 

o 

and fixed frames coincide. All components of torque, 
momentum and inertia are referenced to the axes. 

b. are exterior torques acting on the contents of 
a volume V which is arbitrary. 

c. Quantities h^ and c^^^ depend on the fields of den- 
sity p(x^,t) and relative velocity, u^(x^,t), where p and 
are independent variables subject to constraints of laws of 
conservation and equations of state. 

d. Let y^^ be a coordinate system rotating with the same 
angular velocity, cj, as the x^ system, but with y^ directed 
always along the instantaneous rotation axis. Then, 


0) . 

1 

X . = — y, 

1 u) 3 


(142) 


(D . 


SO that — are the direction cosines of the rotation axis 


0) 


relative to the reference axis; thus, the wobble components are 
0 ), and o)_ . 
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The expression in (141) is based on a rigid Earth. 
From Jeffreys (1970) , the total external gravitational poten- 
tial is 


U = 


Ga' 


+ + 


ka' 


U, 


(143) 


and, from the same source, MacCullagh's formula for the grav- 
itational potential from a distant point is 

T T GM G,- GM ■ G „ ,2^3-5ij. 

U = — + (A + B + C - 31) = — + C. . (r n -3x x-* ) 

'>-,3 r 


(144) 


2r 


2r 


Equating equivalent components of (143) and (144) , 


GC^j (r2ri^^-3x^x^) =2ka^U 


(145) 


i i 

where n is the unit matrix. For a rotating deformable Earth, 
the concern is for the centrifugal potential, 

U = r^ - (u)^x^) (r^n^^-3x^x^ )] . (146) 


Substituting the right-hand part of (146) (discarding the radial 
component, ^^r^, since it has no influence) into (145) yields 

(147) 


^ ka^ 

C . . = U • 0) • , 

1] 3G 1 ] 


which can be written as 

ka ^ 

C. . = In • • + u.o). + Constant , 

ij ij 3G 1 j 

where I = is the inertia of the sphere^ free of rotational 

deformation. This determines the constant and we have 
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+ 


ka^ 

3G 


(o) . OJ . 

1 D 


1 2 ^ 
^ n . . ) 
3 1 ] 


(148) 


Orienting the axis along the rotation axis, a)j=w^ = 0, o)^=a), 
and 

k a^ 2 k a^ 

C-A 

Setting = H, (H is known as the precessional constant) 

we obtain 


k = 
s 


3GHC 

, 5 , ,2 

a oi) 


(150) 


where k^ applies to the rotational case carried on for the 

total of the Earth's rotation time, and is known as the secular 

Love number. It is considered equivalent to the fluid Love 

number, k^. The secular Love number has been evaluated as 

k =0.96 (Munk and MacDonald (I960)), and this value is applied 
s 

to k^. 


From (147) and' (150) , we find the products of inertia 

to be 


c . . = 
13 

k C-A 

5 0) . (JL) . 

2 13 

f CO 


(151) 

from which. 




^13 = 

^ (C-A) m 

Kf 1 

, c - ^ (C-A) m 

2 3 k^ 2 

(152) 


Substituting these into the first two of (140) (for the moment, 
disregarding h^ and N^) , 
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V m 


4 >, = 




m 


In the complex mode, these are written as 


A k , . m. 

<p = = — (m - 1 — ) 
k^ uj' 


We set 

u ^ ^ 


^ I 


( 153 ) 


(154) 


(155) 


Introducing 



o 

r 


(155) into (141), 
” ' *D ■ “-V "" *D 


(156) 


C— A 

The approximation depends upon the magnitude of — and here 
the error 0.1 percent. The expression may be interpreted 
as that part of the excitation function (J) that is due to rota- 
tional deformation. Then, 


im + o^m = 0 


which differs from the corresponding rigid Earth expression 
when <p=0, 

im + a m = 0 , 

r ' 

in that the frequency of free nutation has been reduced from 
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to 


o. 



a 

r 



(157) 


(We note again the adopted value, k=0.29, given in Section 2.9.) 


The excitation function, (|) , may be defined in terms of 
a number of excitation components, dependent upon the nature 

of deformations or other geophysical phenomena. These will vary 
the frequency of wobble, but, if the excitation functions are 
properly defined, the corresponding wobble frequencies can be 
introduced. We consider excitation functions in the form. 


1 



9 


where k is some transfer function separately considered for 
various geophysical excitations. Some excellent examples are 
given by Munk and MacDonald (1960), regarding this matter and 
various rotation poles for igiven excitations . The differential 
equations, using the general format are 

, « 

— + m = S' , (158) 

m^ = , (159) 

for the polar wobble and the Earth's rotation, respectively. 

The solution of (158), which is of the form, 

dM 

— + P(t)M - Q(t) = 


0 
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which can be solved using the integrating factor, e'^^^^=e 
is 

m(t) = - ia^ J'F(t) e"^^o^ dr] , (160) 

0 

where is an arbitrary complex constant. Consider, first, 
'P(t)=JH(t), where H(t) is the Heaviside step function. Then, 
m^=0, m=0 for t<0, and 

m(t) = Jd-e^^flt) ^ t>0 . (161) 

The second case of interest is that for a harmonic excitation 
of any frequency, o: 

t = 4'^cosat + Y^sinat . (162) 


c c s s 

The four numbers T^, 4'j, determine phase and amplitude 

of the excitation. It is also convenient to use alternate 
forms (See Section 2.4). The solution for the harmonic excita 
tion is 


a 4* 

, 10 t , 0 

m ( t ) = m e o + 


Oo-o 


o 4^ 

lot , ^0 -lot 
e + 


0^+0 


(163) 


with the forced frequency, o. 

The solution of (159) is 


m . 


= 4*. 


(164) 
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The effects of polar wobble on observed astronomical 

latitudes and longitudes (or times) are determined as follows; 

The transformation matrix, a.., for the transformation, 

1 j 


IS 


^(Astro) ^ij ^(CIO) 


a . . = 

13 


cos X 

I 

0 

sin X 


0 

1 

0 


-sin X 

E 

0 

cos X 


1 0 0 
0 cos y -sin y 


0 sin y cos y 
P P 


0 

1 


-X 


-Yr 


(165) 


Ax 


where Xp=-^^, linear, positive wobble 

displacements in the prime meridian and 90° E. longitude direc- 
tions, respectively, and R is the Earth's mean radius. The 

small magnitudes of x , y are considered in the approximations 

P P 

employed in (165) . The transformation from the instantaneous 
rotation axis to the CIO axis is given by 


^ (CIO) ^ji ^(Astro) 
or 


COS(|)COSX 


' 1 0 x ' 

P 


cosd) cosA 
^a a 

cos(|)sinA 


0 1 Yp 

= 

cosd) sinA 
^a a 

sincj) 


f^P ”^P ^ 


L 1 


(166) 


( 166 ') 
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Utilizing first-order Taylor expansions. 


sincj) = sincj) + A<|) cos(|) , 
Si 9^ 


coscf) = cosij) - A(j) sincj) 

3. Si 


cosA = cosA - AA sinA , 

Si 3 


where A(j)=(j)-(t) , AA=A-A , 

3 .3 


we eventually find 


({) = (() - (x "cosA + y "sinA )" , 

3 p 3 p 3 

A = A - (x "sinA - y "cosA ) " , 

a p a ^p a 

‘ - yp"=°s^a> "'^“♦a 


(167) 

(168) 
(169) 


where in (167) through (169), x ", y " are given in arc seconds, 

P P 

and the corrective terms are in arc seconds, as well (units are 
seconds of time in (169)). We mention, again, the convention 
difference, geodesy vs. astronomy, with the former talcing y as 
positive toward 90° E. longitude, and the latter, otherwise. 

This difference is reflected in equations given, for example, 
in Bomford (1971) and Mueller (1969) , compared with those taken 
with the geodetic convention, as above. 

The investigation of polar wobble is well documented 
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in the literature; ^ ^ however, despite a plethora of research 
and theorizing, there is no real agreement about the source 
and dissipation of the force which appears to permanently 
maintain the Chandler component of the wobble. Most agree on 
the existence of a rather consistent annual component, and 
some fewer, on an equally consistent semi-annual component. 

The instantaneous polar wobble represents the integra- 
tion over the Earth's lifetime of all mass redistributions 
within, under and above its crust. Periodic contributors are 
seasonal motions of the atmospheric mass; ocean and bodily 
tides, with additional contributions due to asymmetric distri 
but ion of land and ocean areas; ground water, snow and vegeta 
tion changes; ocean effects such as surface height variations 


In addition to the authors already referenced above, 
those reviewed in connection with the EOM program include 
Anderle (1973) , (1975) , Arun and Mueller (1971) , Ben-Menahem 
and Israel (1970) , Beuglass and Anderle (1972) , Capitaine 
(1975) , Chinnery and Wells (1972) , Colombo and Shapiro 
(1968), Currie (1974), Feissel et al. (1972), Gaposchkin 
(1972), Graber (1974) , (1976) , Guinot (1965) , (1972) , Hau- 
brich (1970), Israel et al. (1973), Jeffreys (1916) , (1940) , 
(1968) , (1972) , Jeffreys and Vicente (1957a) , (1957b) , Lam- 
beck (1972), Liu et al. (1974), Mansinha and Smylie (1967), 
(1940), Markowitz (1940), McCarthy (1974), McClure (1976), 
Mueller and Schwarz (1972) , Munk and Grove (1952) ,, Myerson 
(1970), O'Connell and Dziewonski (1976), O'Hora and Thomas 
(1970) , Okazaki and Nasaka (1972) , Pedersen and Rochester 
(1972), Pines and Shaham (1973), Proverbio et al. (1972). 
Rochester (1970), Rochester et al. (1974), Rosenhead 
(1929), Rudnick (1956), Runcorn (1970), Rykhlova (1967), 
(1969) , (1974) , Shimazaki and Takeuchi (1972) , Smylie and 
Mansinha (1968) , (1971) , Smylie et al. (1970), Toomre 
(1974), Vanicek (1969), Vicente and Yumi (1969) , (1970) , 
Walker and Young (1955) , (1956) , Wako (1972) , Yashkov (1965) , 
Yatskiv et al. (1973), Yumi and Wako (1970). 
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(inverted barometer phenomena), although mass transport con- 
tributes only an insignificant amount. Efforts to directly 
analyze these periodic effects fail due to lack of complete 
coverage by observation and processing techniques; however, 
spectral analyses of latitude data, observed over a number of 
decades, afford consistent and apparently quite representative 
excitation functions. Refer to Gaposchkin (1972) and Wilson 
(1975), for example, and to Section 5.2 for a detailed descrip- 
tion. 

Spectral analyses of ILS polar data from 1900 to 1975 
are shown in Figures 2.8 through 2.10. The plot in Figure 2.8 
is unfiltered; the plot in Figure 2.9 represents time data 
treated with a Hann filter; the plots in Figure 2.10 are for 
four groupings of 30-year series against the overall 76 year 
series. This last figure gives the power in logarithmic form 
to provide better low-power detail. It can be seen that all 
plots show a specific high power at precisely the 2tt frequency 
(annual period). Figure 2.8 shows three additional peaks at 
5.21, 5.39 and 5.65 rad/yr, while Figure 2.9 shows a fourth at 
5.03 rad/yr. It might be assumed that the unfiltered indicates 
three Chandler frequencies, and the filtered one, four; how- 
ever, Pedersen and Rochester (1972) suggest that too much con- 
fidence should not be placed on the discrimination validity of 

^ ^ Hann filtering represents a premultiplication of time domain 
data by h ( 1-cos (2iTj/n) ) , where j represents the point number of 
the equispaced time series of n elements. The purpose is to 
cause spurious side bands to decay at an inverse third- rather 
than first-order rate. 





Figure 2.10 Power Spectra (Without Filter) 
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spectral analyses due to the limited (in their case, 70 year) 
observation period, and that the Chandler period may really 
be a single valued quantity. We agree with the statement of 
non-validity, but the matter of whether the Chandler frequency 
is a single or multiple valued function is not clear. The 
many contributors to its value (from oceans, crust and core), 
and the heterogeneous constituency of crust and upper mantle, 
leave this as a highly unresolvable matter. In addition, the 
Q factor (discussed below) has a definite influence on the pe- 
riod, and there may be several Q values involved which could 
have radically different values for oceans and Earth. For 
computation purposes we have arbitrarily chosen a Chandler 
frequency^® of 5.3 rad/yr and a Q factor of 50. (These values 
represent implied consensuses (general averages).) 

Polar plots of BIH data are shown in Figure 2.11 for the 
period 1966.04 to 1968.96. The two plots show the motion with 
and without the annual and semi-annual excitation components. 
The process for removing these periodic components is described 
in Section 5.2. A comparison of the ILS (IPMS) with the BIH 
data (both corrected for annual and semi-annual contributions) 
is given in Figure 2.12 for the same period. 

The fact that the Chandler period is longer than the 
theoretical Euler (rigid Earth) period is due to bodily defor- 
mation and effects of oceans and liquid core. These involve an 

Several authors (including, but not limited to. Chandler 
(1891) , Graber (1076) , Jeffreys (1940) , (1968) , Munk and Mac- 
Donald (1960), Rudnick (1956), and Yashkov (1965)) have de- 
termined one or more Chandler frequencies based on spectral 
and other analyses. 
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Figure 2.12 ILS vs. BIH Polar Plot 



exchange of energy and imply some dissipation rate at the 
Chandler frequency. This requires that the Chandler frequency 
be complex, the dissipation appearing in the imaginary term, or 


d =0(1 
0 0 



(170) 


where the Q, or quality factor, is the inverse of what is 
known as the specific dissipation function, given by 


2jt _ AE 
” E 




E dt 


(171) 


which is defined as the ratio of energy dissipation to total 
energy, occurring in a single cycle. Kaula (1968) shows that 
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2 QL 

the kinetic energy has a phase lag — 




( 172 ) 


and, in terms of thermodynamics, the dissipation factor is 
related directly to the change in entropy per cycle, AS, or 


1 _ TAS 
Q 2ttE 


(173) 


where T is the Kelvin temperature. If, then, in Eq. (160), 

S' (t) were to remain constant, wobble motion would decay in a 

time — and the axes of rotation and excitation would ulti- 
mo 

mately coincide; i.e., noting 

/ (174) 

lim m(t) = lim e^^o^je ^*^m + "F (e ^*^o^-l)] = 'F (e ^*^o^-l) ; 

t-yO t^O 

thus, characterization of Chandler excitation functions re- 
quires frequency and dissipation factors for each component, 
as well as initial magnitudes and orientations. 

Attempts to demonstrate the source of energy driving 
the Chandler wobble have not been too successful. Wilson (1975) 
concluded that meteorological variations account for a signifi- 
cant portion, but not all of the wobble, and that other non- 
meteorological sources, such as earthquakes should not be ex- 
cluded from consideration. Mansinha and Smylie (1967) , (1970) , 
Smylie and Mansinha (1968), Smylie et al. (1970) all indicate 
a high degree of correlation between changes in pole path and 
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earthquakes with magnitude, M>7.5, in recent years. Pines and 
Shaham (1973) claim to show that the Chandler wobble and earth- 
quakes may have a common excitation source with the right sign 
and strength to explain the "pumping" of the Chandler wobble. 
O'Connell and Dziewonski (1976) conclude from analysis of 234 
earthquakes with M>7.8 (occurring between 1901 and 1970) that 
earthquakes represent the major factor in wobble excitation. 
Myerson (1970) finds strong support for the Mansinha and 
Smylie (1967) theory but feels that the earthquakes, themselves, 
are not the source but a parallel effect to the wobble excita- 
tion. Ben-Menahem and Israel (1970) find that an earthquake of 
magnitude 8.5, under favorable circumstances, may suffice to 
maintain the wobble for about one year; hence they deduce that 
earthquakes may at most account for 30 percent of the observed 
secular polar shift. A strong rejection of the Mansinha and 
Smylie theory was made by Haubrich (1970) who found that the 
excitation, from actual historical events, is at least an order 
of magnitude too low to maintain the observed wobble. O'Hora 
and Thomas (1970) report negatively on the earthquake as being 
a significant source, while Chinnery and Wells (1972) feel the 
matter is still open but not resolvable with present data. 
Runcorn (1970) feels that treatment of earthquake effects as 
step functions is incorrect and that the input should be impulse 
types. Israel et al . (1973) feel that eaj-thquakes are insuffi- 

cient to maintain the Chandler wobble. It was determined in the 
course of this EOM study that seismic disturbances leading to 
earthquakes could, indeed, be responsible for the major portion 
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of the Chandler wobble power as well as the secular drift of 
the mean pole. The seismic contribution is impossible to 
realistically preprogram and, since its magnitude far out- 
weighs that of any other possible source (and these issues 
are so thoroughly treated in the literature) , no further dis- 
cussion will be made of Chandler components. Details of the 
EOM findings are given in Section 5.2. 


2.15 


Earth Tides 


The complete tidal potential can be written as 


W = 


GM 


GM 


GM \ 


r(l-2^os^i+— ■ 
r 2 n=0 

r 


(f)" P„(.) 


(175) 


where, referring to Figure 2.5, p=PM, r=x^, R=OP, and P^^ are 
the Legendre polynomials described in Section 2.6. From (53) 
and (175), 


= — , a constant affording no tidal effect, 

GM 

Wj = — Rcosf. , which provides a radial potential applying 
r^ 

to the Keppler-Newtonian equations for orbital motion, but is 
not part of the tidal effect, and. 




GM r2 


(3 cos^C ~ i) 


(176) 


W = ^ (5 cos^f, - 3 cos^) , (177) 

A L 

r ^ 

PM R** 

cos'*;’^ - 30 cos^i; + 3) , (178) 

r ^ 
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Eq, (176) is the same expression given earlier in (75) and (92) 
note, also, that R<<r, which explains the relative vanishing of 
, etc., with respect to . Substituting 


cos<; - sin<()sin6 + cos4>cos5cos (t+A ) 


(179) 


and the decomposition formula,^’ 


P^(cosiJ;) = (cos0 ) (cosO ' ) + 


j ; R (e,t)R (6', A] 

^ v* (n+m) 1 nm , nm ' 
m=i ' 


+ S^^(0,t)S^^(0' ,A)1 ; 

nm nm J 


(180) 


where 4 > / ^ is the latitude and longitude of the point P, respec- 
tively, (S , t is the declination and hour angle of the exterior 
body, respectively, 0=90°-(f>, 0'=9O°-6, 

=nm ' P„„(co30)sln(mt). 


into Eqs. (176) through (178), we eventually obtain. 


W = 


3GMR' 


3 ( sin^ 0- j) (sin^ 6-j) + sin2(J)sin26cos ( t+A ) 


W3 - 


4r ^ 

+ cos^(|)Cos^ (Scos2 (t+A 
GMR^ 


8r 


J , (181) 

2 ( 5sin ^(|)-3sin(t) ) (5sin ^ 6-3sin6 ) + 3coS(()COs6 ( 5sin^ (1>-1) 


(5sin ^ cS-1 ) cos (t+A ) + 30sin(|)sin6cos^4icos^ 6cos2 (t+A ) 
+ 5cos^^cos^6cos3 (t+A)] , 


^ ' Refer to a text discussing spherical harmonics; e.g., 
Ileiskanen and Moritz (1967) . 
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These can be combined into long period, diurnal, semi-diurnal. 


ter-diurnal, quarter-diurnal, etc., tides; 


Vj 16 (3sin^())-l) (3sin^6-l) +^^^{5sin^(l)-3sin(()) (5sin^6-3sin6) 

64r' 


+— (3 5sin‘*(|)-30sin^(|)+3) (35sin*‘6-30sin^6+3) 
2 
r 


24R 


+ |^48sin2(})sin2(S+^;— cos(J)COs6 (5sin 4>-l) (5sin 6-1) 


+10B_cos<J)cos6 (7sin ^4)-3sin(|)) (7sin^6-3sin5)l cos(t+X) 
2 J 


+ [4 8cos^cf)cos^'6+- — ^sin(j)sin6cos^ (l)cos^ 6+=-^^— cos^ (j)Cos^ 6 (7sin^(j)-l) 

^ r^ 


(7sin^6-l)j cos2 (t+X) } 

+ [ -— cos ^ ({icos ^ 6+ ^ sin(j>cos ^ t}>sin6cos ^ 6^ cos3 (t+X) 


cos‘*({)cos‘*6cos4 (t+X) 


(182) 


For all practical purposes of the EOM, the expansion of as 

given in (181) is adequate for W. 

As mentioned in the introduction, the total tide rep- 
resents the sum contributions of all its components. There 
are at least 31 major contributors; each operate like a sep- 
arate tide inducing body traveling at its own individual rate. 
It has been found that all tides can be expressed in terms of 
one or more of a set of six variables. Doodson (1921) de- 
scribes the procedure and sets up a systematic classification 
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for these tidal waves. The defining variables are 
T = 15°. 0 t+h-s-X 

s = 277°. 0248 + 481267°. 8906T + 0°.0020 t2 , 

h = 280°. 1895 + 36000°. 7689T + 0°.0003T^ 

(183) 

p = 334 .3853 + 4069°.0340T - 0°.0103T^ , 

N' = 100°. 8432 + 1934°.1420T - 0°.0021T^ 

Pj = 281°. 2209 + 1°.7192T + 0°.0005T^ 

where^ 

T = (Julian Date^Qj^gj - 2415020 . 5) /36525 , 

X is the local mean lunar time reduced to angle, s is the 
Moon's mean longitude, h is the Sun's mean longitude, p is 
the longitude of Moon's perigee, N' = -N, where N is the 
longitude of the Moon's ascending node, p^ is the longitude 
of the Sun's perigee and X is the longitude of the station. 

The speeds per mean solar day are 

T = 360°. 0 - 12°. 19074939 
s = 13°. 17639673 , 

h = 0°. 98564734 , 

(184) 

p = 0 .11140408 , 

N' = 0°. 05295392 
Pj = 0°. 00004747 

Doodson's classification is in terms of the arguments defined 
in (183) which would apply to a given tidal wave. For example, 
the argument for the Rj wave is 

: 2t - h + 2s - p j 
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Doodson's classification is a six digit one, 

a, (b+5), (c+5), (d+5), (e+5) , (f+5) , (185) 

where the letters, above, are the coefficients of the wave 
expression, 

a + bs + ch + dp + eN* + fpj 
Finally, there is a decimal separating the third from the 
fourth coefficient in (185) . The classification for is, 
therefore, 

R : 274.554 

Note that the larger the argument number, the greater is the 
speed of the component tidal wave. 

The harmonic development of the equilibrium tide is 
discussed in Section 2.16 where its application will have 
more relevance. 

The influence of luni-solar tides on the deformable 
Earth was described in Section 2.8. It was assumed there 
that the tide was an equilibrium one, and no consideration was 
made of what is known as the secondary tidal effect, that due 
to Earth loading by oceanic tides. The equilibrium concept 
applies fairly well to Earth tides with the exception of some 
phase angle in the tidal bulge due to a given tide. Phase 
angles vary (they may be positive as well as negative) at a 
given station for a given tide.^® Although expressions are 
given (see, e.g., Kaula (1968)) for the lag angle, they depend 
upon variables which are not available. 

2 8 


See Melchior (1966) . 




90 


Darwin (1879) appears to have made the first calcula- 
tions on tidal loading of the Earth's surface (secondary 
effect), taking into account the distortion of the Earth's 
surface. The problem of effects on the vertical deflection 
was presented in a rather comprehensive fashion by Lamb (1917) . 
Theoretical effects based on various Earth models were dis- 
cussed by a number of authors including Alsop and Kuo (1964) , 
Caputo (1962), Longman (1962) , (1963) , (1966) , and Takeuchi 
(1950) , the latter being an especially comprehensive treatment. 
Lambert (1974) applied tidal spectroscopy methods (developed 
by Munk and Cartwright (1966) ) to analyses of total gravity 
and tilt data to develop a weighting system to describe the 
Earth tide more accurately in an environment of complex ocean 
tide perturbations. A study by Nishimura (1950) indicated 
that elasticity of the Earth's mantle is 20 to 30 percent 
smaller in the N-S than in the E-W direction. Kuo and Ewing 
(1966) noted that ocean tidal loading effects decay approxi- 
mately exponentially as a function of station distances from 
the effective ocean water; however, for stations in the imme- 
diate ocean vicinity, tidal loading effects on tidal gravity, 
due to land tidal conditions, is very complex, small and 
irregular. Kuo and Jachens (1970) investigated the conti- 
nental tidal gravity profile across the United States and 
found that tidal characteristics of both oceans entered, but 
they could find no observable correlation between tidal gravity 
parameters and regional geology. Lambert (1942) had already 
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determined that at inland stations the secondary effects may 
be perceptible for considerable distances, even as much as 
4000 kilometers, and that at coastal stations, secondary 
effects could mask out entirely the primary effects; his work 
suggested that (See Eq. (181) ) component hour angles should 
be modified by subtracting appropriate phase angles e^, e^. 

The effects of ocean loading can be measured by using 
gravimeters, tiltmeters, horizontal pendulums and extensome- 
ters (strain gauges) . A complete description of instruments 
and techniques is given by Melchior (1966) , and Slichter 
(1970) provides an excellent and comprehensive treatment of 
observational efforts to measure interaction of ocean and 
land tides and describes efforts and results at a number of 
stations. Investigation of local effects have been conducted 
by Beavan (1974), Egedal (1959), Harrison et al. (1963), 
Lambert (1970) and Lennon (1962) . 

Interactions between Earth and ocean tides were studied 
by Kozai (1965) who found that effects of the tides are about 
10 percent of the direct luni-solar effects on the satellite, 
and that eccentricity, as well as the semi-major axis, are not 
disturbed except for short-term periodicities. Lambeck et al. 
(1974) found that consideration of the attraction of ocean 
tides was quite important in utilizing satellites for deter- 
mining the Love number k, and that an equilibrium theory was 
not adequate for correcting for ocean tides; after more appro- 
priate treatment, improved values of k were obtained, but poor 
solutions for phase angles indicated a need for more precise 
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satellite tracking data and tidal models. 

There appears to be no analytical means for introducing 
corrections for secondary effects into the EOM, and no bases 
for adopting values for phase angles. Munk and MacDonald 
(1960) indicate a relationship between phase angle and Q (for 
large Q) , 

e = tan-M^) , 

which, for a Q of 50 (adopted arbitrarily for the EOM) , indi- 
cates an e=0°.57. This, of course, is merely a generalization, 
and values given in Melchior (1966) indicate no possibility for 
presupposing a fixed phase angle. In fact, the wide range of 
magnitudes and signs leaves this parameter as a unique element 
to be specified only from direct observations for given compo- 
nents of the tide. The EOM does not incorporate phase angle 
values for given stations but, instead treats the luni-solar 
influence as an equilibrium tide; however, values as they 
become available, can be readily introduced if desired. 

The number of tides considered for the EOM have been 
reduced to three components in each of the long period, diurnal, 
and semi-diurnal classes. These are the S , M , M^, O , P , 

5a in t , * 

Kj, Nj » M 2 , S 2 components with argument numbers: 

"s : 057.555 0,: 145.555 N, : 245.655 

sa 1 2 

M : 065.455 P,: 163.555 M 2 : 255.555 (186) 

m 1 ^ 

M^ : 075.555 Kj : 165.555 S 2 : 273.555J , 



2.16 Ocean Tides 


93 


The familiar ocean tide appears to an observer on shore 
as a repetitive, systematic rise and fall of the ocean level. 
The range of the tide varies with relative positions of Sun 
and Moon, and, sometimes, quite spectacularly: when the coast- 
al configuration lends itself to resonance conditions,^® which 
can cause extremes in highs and lows. Constraints imposed by 
land systems provide unique local circumstances which can be 
closely characterized by a finite system of descriptive param- 
eters called harmonic tidal constituents. 

The method of harmonic analysis was first developed by 
Thomson (1875) , and perfected by Doodson (1921) into a suc- 
cessful system; Doodson (1957) was also the first to success- 
fully treat for shallow water constituents. In practice, 
tidal constituents are determined by harmonic analysis of ex- 
tended tidal height measurements, in the following general 
fashion; Let the observation equations be written for the k 
constituents as 


k 

y ■ A COSO) t. 

7 ^ ^ 


+ B sinu) t. 
r r 1 





(187) 


where i = -n, -n+1, . . . , 0, . . . n-1, n, are the observa- 
tions taken at uniform intervals At, where n=0 is the middle 
observation of the time series, for which t^=0; A^, are 
magnitudes of the tidal constituents to be determined; hg is 


2 9 


See Garrett (1972) 
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the mean level to be determined; are the known frequencies 
associated with the selected tidal constituents, h. are the 
observed amplitudes at the times t^, and v^ are the residuals 
whose summed squares are to be minimized. As many constituents 
(A,B) ^ are chosen as will completely satisfy the unique char- 
acteristics of the local tide; e.g. , for Anchorage, Alaska, 
Zetler and Cummings (1967) required 114 constituents to prop- 
erly define the tide. 

The height of the tide at any coastal point and time may 
be represented by the expression, 

k 

f = f- + ~'y~' f h cos(o) t-(j) ) , (188) 

0 n n n ^n 

n=i 

where f ^ is the mean height in respect to a reference level; 

h^ are mean amplitudes of constituents during the extended 

period over which h were determined; f are node factors 

depending on variables p, N' and p^ for which the periods are 

approximately 8.6, 18.6 and 21,000 years, respectively; t is 

the time in UT; <b = - qX + k - u S - (V +u ) ; X is the lon- 

n n n n n 

gitude of the point east of Greenwich; q is the species number 

(0 for long-period, 1 for diurnal, 2 for semi-diurnal, etc.); 

(V +u ) are equilibrium constituents at Greenwich; S is the 
n n 

number of hours which standard time at the site follows Green- 
wich; ~ *^n ~ qX-o)^S, and, together with h^, are called the 
harmonic constants for the point. Procedure is to determine 
<})^ for a given period, together with h^ and find the value g^ 
from the above. Averages of h^, g^ over an 18.6 year cycle 
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are utilized, and are phase lags of the tidal constituents 
behind the corresponding equilibrium constituents at Green- 
wich. Constants f and (V +u ) are computed and available in 

n n n 

Shureman (1941) . 

Equation (188) fails for tides in open oceans (pelagic 
tides). Water speed limitations of various tidal waves, and 
Coriolis effect present problems which must be solved in some 
other manner. The basis for modern efforts are solutions of 
Laplace's tidal equations. The involved mathematics and 
restrictions are given by Hendershott (1975) , as well as the 
general status of numerical modeling of ocean tides. Diffi- 
culties involved in modeling the pelagic tide include require- 
ments for specifying ocean depths and boundaries. Since both 
are highly irregular and not analytic, a complete mathematical 
solution is impossible. Investigators must specify boundaries 
in the forms of meridians and parallels, and uniform depths, 
in order to achieve any solution at all.^” 

The resultant of motions of tidal waves affected by the 
Coriolis forces and moderating influences of land masses, 
result in each component of the tide forming individual amphi- 
dromic systems. A sample of an amphidromic system for the 
tide, taken from Luther and Wunsch (1975) is shown in Figure 
2.13. 

Amphidromic points are points where, for that particular 
tidal component, the water level remains fixed in height. 

see Doodson (1927) , (1936) , Goldsbrough (1927) , (1929) , Hen- 
dershott (1975), Pekeris and Accad (1969), Proudman (1916), 
(1925) , (1932) , (1936) , (1944) . 
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Figure 2.13 Amphidromic System 



Radiating out from the amphidromic points are solid cotidal 
lines representing the phase of the tide in even hours from 
the high tide at Greenwich; i.e., the phase lag to high water 
after passage of the representative body over Greenwich (in 
this case, in solar hours) . The dashed lines are cb-range or 
co-amplitude values representing the tidal amplitudes in cen- 
timeters. Computational procedures for amphidromic systems 
are modified by boundary conditions presented by tidal data 
at continental coastal stations, as well as those on island 
systems. These coastal and island tidal data are collected 
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and distributed by the International Tidal Bureau at Monaco. 
(Tidal parameters for about 4,000 stations distributed through- 
out the World are included in the computer package.) 

Additional authors involved in determinations of amphi- 
dromes whose works were reviewed included Cartwright et al. 
(1969), Defant (1961), Fairbairn (1955), Irish et al. (1971), 
Munk et al. (1970), Neumann and Pierson (1966), Parke (1976), 
Platzman (1975) , Zetler et al. (1975) . 

A complete set of amphidromes (for all practical pur- 
poses) would include the six annual and semi-annual compo- 
nents given in (86). The tide for a given point and time 
would be given by the algebraic sum of the reduced levels. 

In practice, phase and amplitude values would be stored on 
tape for specified grid points, and extracted by using double 
interpolations (See Section 3.3.5). 

The long-period tides (S , M , M^) are treated as 
equilibrium tides. From Maximov (1970), the equilibrium equa- 
tions are (with W in centimeters) , 


H(S ) = 0.95 (l-3sin^<j))cos2h , 

Sd. 

H(Mm) = 1 . 08 (l-3sin^(j)) cos (s-p) 
H(M^) = 2.05 (l-3sin^(j))cos2s , 


M189) 


where H is the disturbance in mean-sea level, (f> is the latitude 
of the point and s, h and p are defined in (183). These values, 
when added algebraically to those obtained from the amphidromic 
systems, give the total tide (See Section 3.3.5). 
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2.17 Geomagnetic Field 

2.17.1 General Note; 

The magnetic field existing at the Earth's surface and 
at satellite heights is made up principally of geomagnetic 
components. These include a combination of what is known as the 
dipole and non-dipole field. The dipole and lower order non- 
dipole (below the 7th harmonic degree) are apparently maintained 
by core-mantle interaction, through what has become known as the 
"dynamo effect"; higher-order components may include contribu- 
tions of magnetic materials in the Earth's outer crust. 

Solar storms and solar wind influence the Earth's magnetic 
field, the former, at times, inducing large field changes, but 
such phenomena are not preprogrammable items, and can be evalu- 
ated only from direct observation. Such considerations are 
beyond the scope of this preliminary EOM system and, as a 
consequence, effects of solar radiation have not been considered. 
There is, however, nothing to preclude introduction of such in- 
formation if and when indicated by future studies, since the 
module designs have been kept open-ended. 

2.17.2 Geomagnetic Nomenclature and' Relationships ; 

The conventional magnetic elements are X, Y, Z, D, I, H, 

F; where X, Y And Z are the magnetic field intensity components 
(usually given in gamma units = 10”® Gauss) in the North, East 
and Down directions, respectively; D is the magnetic declina- 
tion or clockwise angle of magnetic north from geographic north; 

I is the magnetic dip or angle below the local horizontal in the 
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magnetic north direction; other elements and relationships are 
F = (X^+Y^+Z^)*^ 

and 

H = FcosI, 

for the total magnetic field intensity and its horizontal compon- 
ent, respectively. The vertical plane through the magnetic force 
vector, X^, is called the local magnetic meridian. It makes the 
angle D with the geographic meridian. Geomagnetic coordinates, 
and their conversions from geographic coordinates are illustra- 
ted in Figure 2.14. 


Figure 2.14 Geographic to Geomagnetic Coordinate Conversion 


N 



B is the magnetic pole with geographic 
coordinates <}>i, Xx; P is the point with 
geographic coordinates (}>, A. Geomag- 
netic coordinates are $, A, and are 
given by 


<J)=sin“^ |^sin(j)sin(|) j+cos<j)COSc}) jCOS ( A-A ^ )] 


The geomagnetic field is conventionally represented^^ 
in spherical harmonics of the Schmidt form (See (49) , (50) , 

(53) for definition). The harmonics have the following inter- 
pretations: 

See Barraclough et al. (1975), Cain et al. (1965) , (1967) , 
Hurwitz et al. (1966), Vestine et al. (1963a) , (1963b) , for 
expressions of the geomagnetic field in spherical harmonics. 
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Axial dipole = Pj , 

Equatorial dipole = pj , 

Centered dipole = pj , 

Eccentric dipole s p^ , pj , , p2 ^ 

Non-dipole field = P^ , Pj » P| ' ^ 3 ' • • • ' 

Not considered = Pg 

The magnetic elements can be developed from the negative 
derivative of the potential or. 


and 


F = -vv = - a .V 

1 


^r = 


iY 

9r 


^0 = 


1 ^ 
r 30 




1 3V 
r sin 0 TX" 


X = - Fn cos 6 - F sin5 
0 r 


Z = F„ sin 6 -F cos 6 , 
0 r 


Y = F, 


S = (j)-ilj = (p - tan ‘ [( 1 -e^ ) tanc))] 


'(190) 


The potential, in spherical harmonics, is 
°° n 

V = a/\(— (g"*cosmA+h"*sinmA) P*'^(cos9) 
r — n r ^ n n n 


n=l ^ m =0 

(a=average radius vector''^6371 . 2km) 


(191) 
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Introducing (190) into (191) , 
n(max) n 


3P"'(0) 


Fq = - V ^(g™cosmX+h^sinmX)-^ 

n=l m=0 


^X = 


F = 
r 


n(max) n 

^ /P^n+2\~' m(g™sinmX-h™cosmX) P™(0) 


2'r> Z. 


n=l 
n (max) 


m=0 

n 


^ (g^cosmX+h^sinmX)P^(6) 


n=l 


m=0 


1^(192) 


where p is given by (15) and, for points at altitude, we use 

(19) for r = p + H, where H is the altitude of the point 

above sea-level. Coefficients g™ and h^ have been selected 

n n 

from Cain et al. (1967) and are stored in the computing package. 
2.17.3 Geomagnetic Field: 

The geomagnetic field appears to have been universally ac- 
cepted as the result of a self-exciting dynamo activated by 
fluid motions of the conducting liquid Earth's core, which act 
like a rotating disc in a disc dynamo. The mechanism is simply 
explained by Takeuchi et al. (1967). The element of interest is 
the westward drift which has been investigated by a number of 
scientists (See Footnote 20 on p. 64) . This element was dis- 
cussed in Section 2.13, and its possible influence on the Earth's 
rotation must be a factor in any further development of the EOM. 

A few comments are in order before closing out this section: 

The dipole is oriented at an angle to the rotation axis of 


These include Bullard et al. (1950), Bullard and Gellman 
(1954) , Inglis (1965) , Parker (1955) . 
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tan [(g!%h;“)Vg jj; an angle of some amount is necessary accord- 
ing to Cowling (1957)^ who showed that a magnetic field symmetric 
about an axis could not be maintained by a symmetric motion. 
Leaton et al. (1965) give the adopted positions of the north 
magnetic dip pole at (|)=75°.5N, A=259°.5E and the south magnetic 
dip pole at <f)=66°.5S, X=139°.4E; they also give the rate of 
change of the dipole moment at (-16+2) r^ gamma/yr, where r is 
the mean radius. Bullard et al. (1950) discuss thermal convec- 
tion in the core and indicate that resultant radial motions 
require that material near the outside of the core rotate with 
a lesser angular velocity than that inside, in order that angu- 
lar momentum be conserved. Yukatake (1972) suggests that the 
angular velocity of the mantle changes with conductivity var- 
iance of the lower mantle (which is important if a small amount 
of leakage of the toroidal field occurs from core into mantle) ; 
he further notes that a one percent change in the dipole moment 
may cause a variation in Earth's rotation rate of 5xl0~^^ sec"\- 
thus leaving open the possibility that the dipole could be a 
cause of short-period fluctuations in Earth's rotation, par- 
ticularly so if toroidal field leakage participates in the coup- 
ling. Rochester (1960) considers the moments of inertia of core 
and mantle and determines that a decrease in the rate of west- 
ward drift by one-tenth of its steady-state value lengthens the 
day by 1.3 msec. Equations relating the westward drift have 
been suggested by James (1968), Nagata (1962) , (1965) , Rikitake 
(1966) and Winch (1968) . 
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2.18 Gravity 

Use of Earth satellites and improved observational accu- 
racies have afforded, not only means, but reasons, for more 
extensive knowledge of the Earth's gravity field. The gravity 
field at altitude is, effectively, an upward continuation of 
the Earth's surface gravity. Differences between observed 
gravity and normal gravity (See Section 2.7), taken at ground 
level, are high-frequency phenomena and can be characterized 
only in discrete form. The situation becomes more analytic 
with altitude, and methods of spherical harmonics become quite 
effective and useful there in describing the gravity field. 
Zonal harmonic coefficients can be obtained from long-term 
orbital observations by analyzing perturbations and secular 
changes in orbital elements; tesseral coefficients can be 
determined only with precise tracking of orbital oscillations. 
The keys to comprehensive results are uniform distribution of 
observational data, a variety of satellites with appropriate 
orbital inclinations, provisions for properly treating atmo- 
spheric drag, and accurate and varied observational systems. 

A great amount of ground gravity data has been accumu- 
lated and utilized in determining the geopotential configura- 
tion and the parameters of a best-fitting representative 
ellipsoid. Data holdings (1959) are given by Heiskanen and 
Moritz (1967) , but the collection has been continuing, subject 
only to physical and political constraints. Procedures for 
utilizing ground gravimetric data for geoid determinations are 
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described extensively in the literature . ^ ^ These data have been 
advantageously combined with satellite gravimetry, providing a 
complementary filtering effect. Results have been reported by 
a number of authors.^** 

Satellite altimetry offers a means for directly deter- 
mining the ocean geoid, and thus the gravity field. Mourad 
(1975) reports on pending experiments, as well as the first 
effort on Skylab in 1973 (McGoogan et al. (1974)). Additional 
reports indicate considerable promise for this method. The 
pelagic tides will enter into this approach in the following 
manner; theoretically, satellite altimetry (SEASAT) could pro- 
vide means for completing the amphidromic characterizations, 
but the tidal influence on the satellite path could render this 
somewhat of a "bootstrap” operation, unless the SEASAT could be 
monitored to isolate these perturbations.^® 

Other satellite gravimetric methods are mentioned by 
Chovitz (1975) , involving satellite-to-satellite tracking tech- 
niques (Comfort (1973) and Yionoulis and Piscane (1972)). A 
satellite application of gravity gradiometry is described by 
Forward (1973), and it is compared with other systems by Glaser 

See, e.g., Bomford (1971), Caputo (1967), Heiskanen and 
Moritz (1967) , Heiskanen and Vening Meinesz (1958) , Stokes 
(1849) , Uotilla (1960) , Zhongolovich (1957) . 

See Caputo (1967), Cook (1963), Gaposchkin (1974) , (1975) , 
Hopkins (1972), Koch (1974), Lerch et al. (1974), Rapp (1971). 

These include Brown and Furry (1973), Lundquist and Gia- 
caglia (1972) , McCandless (1974) , Weiffenbach (1972) . 

^® See Musen and Estes (1972), Musen and Felsentreger (1973). 
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and Sherry (1972) . 

Limitations of spherical harmonic determination and va- 
lidity of expression are discussed by Gaposchkin (1975) . He 
suggests that satellite data are strong only for degrees less 
than 12, and spherical harmonics are probably reliable only 
for degrees less than 10. This is in general keeping with an 
empirical relationship given by Kaula (1968) : the magnitude of 
normalized coefficients of the Earth's gravitational field up 
to about degree £=15 is approximately Despite these 

limitations, it was decided that the gravity field of Gaposchkin 
(1974) , in spherical harmonics to the 18th degree and order, 
should be utilized to those extents in the EOM. 

Spherical harmonics for gravity have traditionally been 
the fully normalized ones shown in Eqs. (51) through (53) . 

Gravity is given (for Earth-bound points) by 
18 n 

(n+1) (^)'' 
n=0 m=0 

where; 

p = jja^cos^6+b^sin^3]^ = a [cos^ 3+ (1-e^ ) sin^ 3^ , 

3 = tan~‘ |^(l-e^ ) *^tan<J)j + 1" .71Hsin2cj)xlO“'* , 

0 = 90°-<J)', r = p+H, = (|)+1" . 7lHsin2(j)xlO“'* i 

a and b are the Earth's semi-major and -minor axes, respectively 
(a = 6378160 meters, b = 6356774.719 meters), GM is the product 

of the universal gravitational constant and the Earth's mass 

/ ' 


B C^'cos mX +S™ sinmxlp'"* - ^^r (1-P //5) , (193) 
n n J n 3 ^ 
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(GM = 3.986010 x 10^“* m^sec"^), (p,BfX) are the ellipsoidal co- 
ordinates of the point, 4) and X are the latitude and longitude 
of the point, respectively, C™ and S™ are normalized harmonic 
coefficients, H is the altitude of the point in meters, P'™ are 

fully normalized Legendre polynomials, to is the Earth's rotation 
rate (to = 7.292105 x 10“®sec“*). For satellite points, the 
right-hand term in (193) is disregarded, and GM = 3.986013 x 10^ '* 
m^sec"^ to include the atmosphere. 



Chapter 3 


EARTH AND OCEAN MODULES 
3.1 General Design 

Each modular component of the EOM has the facility for 
ready update when it is desired to change or increment the 
system. The main program is known as the Input-Output pro- 
gram (lOP) , and is the place where the user selects a par- 
ticular operational mode, introduces data to be processed, 
and necessary input-output coordinate transformations are 
made. Details of mode selection are given in Section 3.2. 

All constants and other parameters which apply to indi 
vidual geophysical modules are put into Common statements 
where possible; data changes are made in the modules where 
they are introduced. A description is made in each of the 
sections giving geophysical module details (Sections 3.3.1 
through 3.3.8) of each parameter identification code and its 
location. 

In principle, the lOP would have a direct tie to each 
geophysical module (GM) as shown in the block diagram in 
Figure 3.1, and in the flow diagrams in Figure 3.2. It was 
decided, however, for convenience, to tie only four of the 
modules to an lOP, and to maintain the remaining four'" as 
individual units. Factors bearing on this decision are dis- 
cussed in Sections 3.2 and 3.3. 
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Figure 3.1 Block Diagram of EOM 
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In Figure 3.1, the (X^) ^ input to the precession module 
defines the epoch of choice for the inertial reference system; 
inter-module data flow' is shown by two-way and one-way symbol- 
ism. The first two diagrams in Figure 3.2 show the input and 
output stages, respectively; interfacing points and are 
keyed to all GM programs. Direct data exchange, not involving 
the lOP, are for interface points through respectively. 

Inter-module operations requiring iterative interactions are 
shown in the third and fourth diagrams in Figure 3.2. 


3.2 Input-Output Modes 

The flow diagram of the system of four modules tied into 
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Figure 3.3 


Flow Diagram of Four-Module System 
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the lOM is shown in Figure 3.3. The system is designed to 
accept station coordinates in either rectangular geocentric 
(metric units) or geodetic (sexagesimal format) with eleva- 
tions in meters. Output, for all cases, is in the geodetic 
format. Additional options can be readily introduced, but 
it is felt that for such isolated cases as might occur, the 
simplest procedure would be to transform them in a separate 
program. The probability of using the geomagnetic field and 
gravity modules exclusively for certain operations led to the 
decision to keep these separate and intact; they are not 
directly related to the lOM, but the lOM can be used to con- 
vert rectangular coordinates to the required geodetic values 
for the modules. The conversion system built into the lOM 
can be reproduced for either, or both, of these latter two 
modules . 

The precession module has little relationship with 
other modules, and its complexity and time-consuming compu- 
tational characteristics indicates a separate role as the 
best practice. 

The polar wobble module requires an elaborate update 
procedure which makes its integration with other modules an 
impracticable concept. Specifics are discussed in detail in 
Sections 3.3.2 and 5.2. 

Input data and updating procedures are discussed in 
the individual GM sections. 
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3.3 Geophysical Module Details 
3.3.1 Precession Module: 

The precession module (EOMPMA) includes 11 subroutines 
in the calling sequence indicated in Figure 3.4. Data and 


Figure 3.4 EOMPMA Subroutine Calling Sequence 
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parameters are located as follows: 

EOMPMA Function is Main; no data. 

INITIAL Accepts read-in data. These include: 

NPRT = Number of print intervals; 

lEPOCH = Epoch of astronomical constants. It is either 
1900 or 1950 (Default) . (In Format 502 card, a non-zero element 
in column 10 calls for 1900.) 

TSTART = Julian starting date. 

DELP = Print intervals. 

TP = Print interval cutoffs. 

AO(l) = Initial Theta value (in radians) for TSTART. 

AO(2) = Initial Psi value (in radians) for TSTART. 

AO(3) = (Leave blank). 

XMO(I) = (Leave blank) . 

CM = Polar moment of inertia. 

AM = Equatorial moment of inertia. 

RR = Adopted constant for Earth's rotation. 

AA = Earth's semi-major axis. 

GK = Gravitational constant. 

XMASS(l) = Sun's mass. 

XMASS(2) = Moon's mass. 

(These cards are punched and located in data section. All 
above data preceding CM must be updated for user's option. 

Item CM and those following may be kept at the user's option.) 

BLOCK DATA Not a subroutine. Contains data not subject to 
change . 

EPHEM Parameter NTTP = 366 indicates maximum number of days 
permissible for a run. Routine reads epheraeris data. These 
data are those provided by the US Naval Observatory converted 
into EOM format. (Coverage: 1970 - 1980.) 

SOLVE Coordinating program for data interpolations, integra- 
tions, and associated computations. No variable data. 

RKVSl Integrating routine. No variable data. 

RKM , LOOKUP , HEAVEN , CURVE Concerned with interpolation and 
application of ephemeris data. No variable data. 

PHI , M1M2 These are used to calculate the Euler angle Phi 
and the direction cosines, respectively. Direction cosines 
indicate diurnal motion of the pole. No variable data. 

RITE Controls the output printout. No variable data. 
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Program outputs are the Euler angles relating the 
ecliptic of 1950.0 to the Earth's axis of figure, and the 
diurnal luni-solar axis motion in arc seconds. 


3.3.2 Polar Wobble Module: 

The polar wobble module is a single program which 
develops a simulated model of the polar coordinates at 
mid-monthly points. Data are output in punched card form 
as well as in print. Simulation is a curve fitting tech- 


nique using the equation, 
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where all quantities are complex except for t^^; mj^ are the 
polar coordinates; m^, are the initial coordinates at time t^; 
4'^ are excitation functions (with subscript a and sa referring 
to annual and semi-annual, respectively); o is the adopted 
value of the Chandler frequency (taken here as 5. 3 (1+0. Old) ) ; 

S is the secular drift in the pole; are Heaviside step 
functions introduced at times, Xj^, and m is the number of step 
functions. Step functions and other elements of Eq. (194) are 
obtained from supplementary programs and operations as follows 
Polar motion data of the IPMS and ILS are not provided 
at convenient equal time intervals which are necessary for 
Fourier analysis. These are converted using WBBLCONV, an 
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adaptation of a cubic spline curve fit. Output is in punched 
cards as well as in print. 

Transformation from the time to frequency domain and 
return is through use of a Fast Fourier Transform (FFT) pro- 
gram called FOURT. Utilizing FOURT as a subroutine to programs 
POWER, REVTRFM AND SPECTRUM, we obtain the following: 

POWER provides the power spectrum of the polar wobble 
time series. The annual and semi-annual power contributions 
are noted and reduced to the average levels in those frequency 
regions . 

REVTRFM returns the frequencies with adjusted powers to 
the appropriate time series. This provides polar wobble data 
with annual and semi-annual effects removed. Both POWER and 
REVTRFM outputs are in punched cards as well as in print. 

SPECTRUM is, in general, the same as POWER with some 
modifications including Hanning and a logarithmic power out- 
put. Output is in punched cards as well as in print. 

ATMOSPH calculates the annual and semi-annual excita- 
tion functions from the REVTRFM output, for various Chandler 
frequencies and dissipation function values. It also calcu- 
lates the annual effects in terms of X and Y components. 

Output is in print. 

/ 

SUBATMO removes the determined annual and semi-annual 
excitations from any polar wobble time series. Output is in 
punched cards as well as in print. 

CWBLANALY and CWBLANALY2 are used to determine the 
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magnitudes and directions of the step functions required to 
duplicate the actual polar motion. CWBANALY is used for the 
initial investigation; when a number of step functions have 
been determined, further operations must be through CWBANALY2. 
The procedure is to determine the points in time which mark 
the limits of a series of points which are generally equi- 
distant along a smooth arc, for each particular section. 

These times are introduced into the program, and corresponding 
complex step values are determined. (Note: When continuing 

with a subsequent series, the last two or so previously 
determined step values should be discarded. The action to 
take will be quite obvious from the punched or printed output.) 

CORRPLOT is used to determine the mean position of the 
pole for given periods over the extent of the polar motion 
time series. Output is in print. 

3.3.3 Earth Rotation Module: 

The Earth rotation module is the first option of the 
lOM routine. The input procedure is the same for all four 
and is as follows: 

The name of the multi-program is EOMIO. The first 
data card is in Format 502 and takes the values: 

NLC = Number of stations. 

NTI = Number of intervals. 

DELP = Spacing between intervals (time) 

TSTART = Starting time (Julian date with a decimal 
equal to .5; i.e., a UT = 0 hours.) 


I 
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PEPOCH = Epoch for the plate tectonics calculation given 
in Julian date. 

BIHIO = Tabular value of UT2-TAI taken from BIH smoothed 
data (pink sheets in the annual report) for the TSTART value. 
When only the Earth's rotation data is required, the value of 
DEEP should be taken as 5.0. Note, also, that TSTART should 
be a multiple of 5 for most convenience in using this section 
in an individual mode. 

SECROT = Secular variation of UT2-TAI determined from 
the most recent smoothed BIH estended series. 

MNU(I) are default options for introducing or excluding 
one or more of the lOM modules. A non-zero quantity in col- 
umns 73,74,75,76 introduces EOMERA, EOMETA, EOMOTA, EOMPTA, 
respectively, into operation. 

The second and subsequent (1 through NLC) cards are 
station description cards. The first element (column 1) 
defines the type of coordinates. A "1" indicates rectangular 
geocentric coordinates. Any other option indicates geodetic 
(sexagesimal) coordinates. After the option is noted, the 
card is read in its entirety in an appropriate program section. 
Data read, in turn, are: 

ICO = Coordinate option. 

LOC = Station name. 

(Coordinates) = Either X,Y,Z (option "1") or latitude, 
longitude and elevation (default option). (H is the elevation.) 

ITPN = Tectonic plate number. 

The Earth rotation module parameters and constants are 
considered firm and should be corrected only as provided in 
the EOMIO input. 


3.3.4 Earth Tide Module: 

The Earth tide module includes 6 subroutines in the cal- 
ling sequence indicated in Figure 3.5.’ Functions are similar 
to those subroutines of EOMPMA of like name. HEAVER is anal- 
ogous to HEAVEN, etc., with the function of RITE being assigned 
to EOMETA. Constants are fairly well established within these. 



120 


Figvore 3.5 BOMETA Sxjbroutine Calling Sequence 



except for XLH = 0.59, XK = 0.29, XLL = 0.07, which are the 
Love numbers h, k, I, respectively, which the user may care 
to vary. 

3.3.5 Ocean Tide Module 

The ocean tide module is, of necessity, incomplete 
until amphidromic data for all six principal (semi-diurnal 
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and diurnal) tides become available. Presently we have in- 
corporated a preliminary set of Mj amphidromic data to set 
the procedure for future availability of these data. The 
procedure for expanding the system is quite obvious from the 
format. The procedure for defining data points is through 
use of a rectangular mercator system of 6° grid intervals. 

The system has a Northing of 102°, a Southing of -120°, an 
Easting of 282° and a Westing of -72°. There will be a total 
of 2280 points for which phase angles and amplitudes are 
given. Conversion of geodetic coordinates to this system 
are given by 

N|^ = ln(cot(45° - , (195) 

E = X - 180° , (196) 

m 

where 4> , X are the geodetic coordinates with the additional 

provision that X is taken as positive from 0° to 360°. 

Constants are fairly firm with the exception of values 

VI (I), which are amplitudes for the S^, S , M and long 

3. S3 in r 

period tides (in mm) , which may be varied at the discretion 
of the user. 

3.3.6 Plate Tectonics Module: 

The plate tectonics module constants for plate motions 
are not considered firm, but they are fairly representative. 
These rotation parameters are given as W(I,J), in units of 
10~^ deg/yr, and are identified by number in Table 2.4. 
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3.3.7 Geomagnetic Module: 

The geomagnetic field module EOMMFA treats the field 
in terms of spherical harmonics using coefficients developed 
by Cain et al. (1967). The input data are: 

LOG = Station name. 

PHI = Latitude (degrees and decimal.) 

XLAM = Longitude (degrees and decimal.) 

H = Elevation in meters. 

MAXN = 10 = Maximum degree of harmonics used. This may 
be taken as less if desired. 

LFN = Type of Legendre polynomial. For EOMMFA this is 
type 3 (Schmidt partially normalized form) . 

IHAR = Option for ellipsoidal harmonics. Leave blank. 

NEWAB = Default option. First station card carries a 
"1"; is followed by the harmonic coefficients, and these are 
followed by remaining station data cards with NEWAB blank. 

Outputs are the magnetic field parameters described in Sec- 
tion 2.17.2. 


3.3.8 Gravity Module: 

The gravity module uses spherical harmonics, with coef- 
ficients developed by Gaposchkin (1974) . Input data for 
EOMGMA is the same as for EOMMFA except that MAXN = 18, 

IHAR = 4 (fully normalized form), and there is no default 
option. Outputs are the gravity at the indicated altitude 
in two forms: G, which is the gravity of a body rotating 

with the Earth, and, GP, which is the gravity for a body 
independent of the Earth's rotation. 



Chapter 4 


APPLICATIONS OF THE EARTH AND OCEAN MODEL 

4.1 General Note on Applications 

There are several areas where in-depth studies of geo- 
physical phenomena can be advantageously conducted using the 
EOM format; however, the complexity of Earth and ocean dynam- 
ics preclude, or at least, inhibit, an effective predictive 
application. Data measurements made with already achieved 
high precisions provide a wealth of material for evaluation 
of cause and effect. Areas of best application are Earth's 
rotation and polar motion, the geomagnetic and gravity fields, 
and the pelagic tides. 

Effective time services provided by the BIH and improved 
timing and rotation monitoring systems can be correlated with 
atmospheric and other disturbing sources to develop more pre- 
cise estimates of secular and aperiodic time changes. SEASAT 
provides an excellent opportunity to establish the pelagic 
tidal components in an empirical rather than the presently 
employed, theoretical, but unnatural procedure. 

There appears to be a less than coincidental relationship 
between the drift of the geomagnetic field and variation of the 
Earth's rotation rate. Correlations could be detected and ex- 
ploited through an EOM application. 

Polar motion studies can be effectively made using the 
EOM and some work was done in this area while conducting this 
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study. The purpose was to test an hypothesis advanced by 
Smylie and Mansinha (1968), and others, that the irregular 
and large variations of the polar path from its theoretically 
proper trace were due to readjustments of Earth's mass associ- 
ated with large earthquakes. The report of this investiga- 
tion is given in the following section. 

4.2 Polar Wobble Study 

Details of polar wobble have been given somewhat exten- 
sively in Section 2.14 and will not be repeated here. The 
polar motion is mathematically described in Eq. (194) , and, 
given enough Heaviside step functions, any polar path can 
be faithfully reproduced. The problem lies in assigning 
causes for these steps, to the amplitudes and frequencies (of 
occurrence) necessary for the fitting configuration. 

The atmospheric influence was rather thoroughly tested 
by Wilson (1975) and considered inadequate as the motive 
force. A number of authors have utilized Fourier techniques 
to determine the atmospheric influence. The consensus was 
that it was a highly predictable phenomenom with a remarkably 
constant annual (and to a lesser degree, semi-annual) period, 
but not sufficient to cause the observed effects. 

An alternative has always been to assume errors in the 
astronomic observations and reductions; however, the ILS 
observations conducted under rather stringent conditions over 
a period exceeding 75 years, have exhibited internal consist- 
encies not warranting the noted excursions. The ILS procedure 
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to utilize the same stars at all stations for their observing 
program precluded significance of systematic star position 
errors. Some systematic error problems still existed, pri- 
mcirily that of evaluating instrumental constants. In the last 
two decades, the BIH has made a significant contribution to 
polar motion monitoring. More recently, doppler observations 
on artificial satellites (by the DPMS) have provided polar 
motion data, and these are combined with the BIH in their 
reports. Although ILS data is not completely consistent with 
BIH data, they exhibit the same trends, and the indication is 
that some factor other than observation error is involved. 

The only remaining "culprit” would be some fairly sudden 
changes in distribution of Earth's mass; possibly signalled 
by a large earthquake. This logical premise has been investi- 
gated by geophysists and seismologists with varying results. 

The problem hinges on finding enough mass shift to cause the 
looked-for results. The consensus is that phenomena associated 
with earthquakes are insufficient by a factor of about 10 to 
account for the polar wobble transients. 

In connection with the EOM development, a test was made 
to determine whether an approximation to the polar wobble 
could be developed assuming that influences took effect only 
at times of very large earthquakes. Earthquakes and their 
dates were taken from O'Connell and Dziewonski (1976) and 
involved 30 events from 1901 through 1964. Magnitudes were 
8.4 or larger on the Richter scale. Using a program called 
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WBLSMOOTHXY, at earthquake times, the Heaviside step values 
were determined to best fit ILS data (corrected for atmospheric 
effect) by a smooth spiral (exponentially decaying due to a 
dissipation factor of 50) joining fixed points. Values of the 
determined Heaviside step values are shown with other data in 
Table 4.1. The program output included a new set of polar 
data points which, when plotted, were reasonably consistent 
with plotted ILS values 

The average earthquake effects (shifts in the center of 
rotation of polar axis from one earthquake to the next) are 
11.9 in X and 12.9 in Y (Units are 10 ® radians) or a diagonal 
value of 17.5, or 36 msecs of arc. Values determined by 
O'Connell and Dziewonski averaged to 19 msecs of arc. The 
determination was within a factor of 2 of the seismologists' 
estimate, well within the uncertainty of their deductions. 
There was no correlation between direction of step values, 
but their estimates were developed from rather uncertain pro- 
jections to subsurface faultings. Comparisons of ILS and 
the simulated earthquake traces are shown in Figures 4.1 
through 4.3 and are representative of the group. The units 
in the figures are 10“ ® radians. 

The conclusion which might be drawn from this experi- 
ment is that it is quite possible that the source, inducing 
variations in radii of the polar motion trace, may be 
mass redistributions in the Earth signalled by very large 
earthquakes. 



127 


, , .artnquaKe maucea Pola. Wob.le Steps 

Table 4.1 Earrnq 


EatthquaXe Interval 
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TO 


1901.71 
1902.79 
1905.29 
1906.12 
1906.62 
1910.46 
1911-46 
1914.96 
1917.46 
1919.38 

1923.04 
1929.21 

1934.04 
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1940.04 
1940.46 
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1941.96 
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1964.29 


Mean Polar 

Coordinates 

X 


Earthquake Effects 
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Figure 4.1 Polar Wobble Plot 1 
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Figure 4.3 Polar Wobble Plot 3 
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Chapter 5 


SUMMARY OF THE EARTH AND OCEAN MODELING PROGRAM 

The EOM in its present state is a first-order approach 
to a geophysical model of the Earth's dynamics. Period of 
effort was inadequate to effectively exploit interactions 
between, for example, the ocean and Earth tides; i.e., the 
secondary effect of the ocean tide. The text of the report 
provides a fairly concise, yet reasonably complete overview, 
of the problems involved in Earth and ocean modeling. The 
programs should be checked more completely than was possible 
in the allotted contract time. It should certainly be con- 
tinued and adapted to some active experimental roles. 

Programs and data are on magnetic tape. In addition, 
listings, with sample data inputs, are provided for the con- 
venience of the user. Approximately 4,000 tidal stations 
(data from the International Hydrographic Office at Monaco) 
are in punched card form. These could be applied as boundary 
conditions in the event that it was decided to investigate 
amphidromic systems. 
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